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I. DESCRIPTION OF MAIN PROGRAM 


A. Program Name 

The name of the main program is TRIPOT, meaning a pot ential prob- 
lem solved by a non-uniform tria ngular mesh. 

B . Functions of TRIPOT 

1 . Generation of a Non-Uniform Triangular Mesh 

The non-uniform triangular mesh is generated on computer by 
the Winslow method (Ref. 1). The essence of the method is to use a trans- 
form to a logical plane, r| which is related to the physical plane x, y by 


9li + 8ii =0 , iiiL + ihuo 

9x 2 Sy 2 9x 2 9y 2 


(la.b) 


and by interchanging dependent and independent variables 


ax + Bx t + Yx tc = 0 

m * 5*1 U 

(2a) 

a yrm + + w# = 0 

(2b) 


where 


a = x? + 
§ 

i 

(3a) 

P = x p x e 


(3b) 

7 = x 2 + 


(3c) 


On the logical plane r|; the domain is composed of equilateral triangles only. 
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Equations ( 2 a, b) are solved by specifying boundary values of 
x, y on the selected logical domain £, rp For fuel sloshing in an axisymmet- 
ric tank we chose a simple parallelogram as our logical domain (Fig. 3 , 
Appendix n). * The interface of arc length s^ is located on i = 1 , j = 1 to N. 
The centerline of length S2 is between the interface and the tank bottom. 

The wall is divided into two parts; the ’'bottom" of arc length S3 and the 
"side" of arc length S4. For cylindrical tanks, the side and bottom are dis- 
tinct; hence it will be used as such and not be redefined, and then and s^ 
are divided into M parts, while S3 is divided into N parts as s-^. For a tank 
such as a spherical tank, the wall is arbitrarily divided into S3 and s^ 
approximately in the ratio s-^ to s^. Parts of approximately equal length can 
be specified. For fuel sloshing, especially with a "folding" interface, more 
accuracy near the contact point and the folding part may be desirable (as the 
machine time and the storage locations increase rapidly with number of 
points on the interface and the side wall). Then one may use a coarser net 
near the center on the interface and away from the contact point on the side 
wall, while s^ and S3 are divided m parallel with S4 and sj, respectively. 
One may also use a finer net near the junction of S3 and s^. In general, it 
is preferable to generate nearly equilateral triangles and to avoid obtuse 
triangles (not necessarily eliminated in the Winslow method). For a folding 
interface, relatively high accuracy of the interface points is needed, and can 
be fed m from outputs of the interface program, SSHAPE, meaning surface 


^All figures in the theory are given in Appendix II. 



shape (see Appendix I). For example, one may use N = 1 7 in SSHAPE (there 
are 16 intervals but take 5 double intervals near the center and 6 single 
intervals elsewhere), dhenin TRIPOT use N = 12 (there are 11 intervals). 
Likewise, a net of N = 24 may be obtained with piecewise uniform spacing. 

With boundary values of x, y specified, the interior net points 
are first solved by a ’'linearized" approximation (Ref. 1), namely 


+ 0 , + 

0£ 2 9p 2 0g 2 0q 2 


(4a, b) 


Then the technique of combined under relaxation of coefficients and overrelax 
ation of dependent variables is used to solve Eq. (2a, b) starting with solu- 
tions of Eq. (4a, b). For Eq. (4a, b), an initial relaxation factor of WA =1.8 
is used. For Eqs. (2a, b), an initial overrelaxation factor of WB = 1. 0 is 
used. The over relaxation factors are improved as given in reference 1. 

For example. 


7& r . , = ETAX = 

'Winslow 


XXf„n + 1 _ n \2 

1 r*ij kT 

II(xP. - x? 1 . " 1 f 
i j iJ 


1/2 


^Winslow LLX 


G0 n 4 T) n - 1 


Winslow 


(5) 


WX 4 ETAX - 1 
WX * n/ ETAX 


(6)T 


tWX is the overrelaxation factor. 
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(“opt) Winslow 


= WAX = 


TT72 " “o 

1 + (1 - X 2 ) 


1 + n/1 - ELX#*2 


- 1 Winslow 
- WO; WO = 0. 01 


/ n + 1 \ 

7 Winslow 


= WX = [j3co< 3pt + (1 - p)o> n l w , 


Winslow 


= RHO * WAX + (1 - RHO) WX; RHO = 0. 05 


( 7 ) 


( 8 ) 


It is noted that if ETAX > 1, use WAX = WX (previous) instead of Eq. (7). 
Let 


XTEMP = [ C 1 X P_ 1( j + C 2 xf_ hj . j+ C 3 x? j _ j 

+ C 4 xf + j, j + C 5 xJ + i, j + l + Cfe*?, j + ll /S'™ 0 (9) 

6 

SUMC = I C k (10 ) 

k = 1 

RX = [xf - XTEMP] * WX (11) 


then 


x. 


n + 1 _ 

i» j 


n 


= x- 


l *J 


- RX 


( 12 ) 


It is noted that in the present problem, we use a simple parallogram in the 
logical plane (Fig, 3, Appendix II). Two examples of constructed nets in a 
physical plane are shown in Figures 4a and 4b, Appendix II. 

2. Sequence of Calculations 

The theory using auxiliary characteristic functions is given 
dijj- 0»jj 

in Appendix II. They satisfy — = XaP on F, — = 0 on W and V 2 ^ = 0 in Vt . 

3n J 3n 
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In our program, vjjj is found by the influence coefficient method. The influence 
coefficients are Fjk and, therefore, one has 

b = F jk if = on F (13) 

or 

<[F] - ~ [i] ) { 4 *} = 0 on F (14) 

which is a standard eigenvalue problem; the k^ 1 eigenvector being ^k* an< ^ the 
corresponding Xk are printed as eigenvalues. The column of Fjk is con- 
structed with a unit k^ element of and zero normal derivative elsewhere. 

3n 

The value of k varies from 2 to N while +(1> 1) = 0 and - = 0, at the 

3n 

center of F. i{j is governed by 

V 2 i|> = 0 (15) 

A finite difference formula of Eq. (15) on the triangular mesh was given m 
reference 1. Let 

A ij = the area of ijtk dodecagon inside the fluid domain 
r. = radius of the ijth point 

w k = 2 (X k 7 k cot d k + \ - l Y k - 1 cot °k>» k = 1 to 6 
(see Fig. la, b)t 

k 3 ' 1 j k k + 1 ' 

t6k> *k can be expressed in terms of tk» Uk> tk - 1 , Uk - 1 and sk» 
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For an interior point {i, j), 


b 2 
£ w k ^k- 40-— ^4 = 0, m=l 


(16) 


with 


He = 1 


For an interface point ( 1 , j), 

6 


I w k M. k - +> - Ay* + (|i)^ [I Sj + i S6 ] r Xj = 0 


( 17 ) 


k = 1 


with 

X^ = X^ — X^ = 0 , X^ — X^ — X 5 — 1 , i - 1 
til 

For ij point on the tank wall, 

6 2 

X w k (^ k “ +) ~ — = °. m = 1 ( 18 ) 

k = l ij 

with 

X3 = X.^ = Xg = 0 and Xj = X 2 = X^ = 1 on bottom wall, 
i. e. , i = M, j = 1 to N 

X 4 = X5 = X& = 0 and X^ = X^ = X^ = 1 on side wall, i. e. , 

/ 

j = N, i = 1 to M 

Xj = X 2 = 1 and X3 = X 4 = X5 = X^ = 0 at corner point, 
i = M, j = N 
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On centerline, 

4* = 0 for m = 1 (19) 

At contact point, i = 1, j = N 

» k (+ k -+>~V + (!±) ({s 3 ) rij = 0 (20) 

k = 1 J ’ 

with 

A.g = 1 and X.4 = X.j = A. 2 = A.g = A.g — 0 

Equations (16) to (20) are solved by the over relaxation method through use 
of subroutine SOREL (1), which contains the same procedure and factors as 
those used in solving for x E see Eqs. (8), (12)] . 

After the influence coefficient matrix is determined, Eq. (14) 
determines the eigenvalues and eigenvectors (the auxiliary characteristic 
functions on the interface). The values of on the wall are obtained by 
the S. O. R. method through subroutine SOREL (2). 

Then subroutine SOREL (1) is again used to solve the boundary 
value problem of a capped fluid (rigid interface) under pitching oscillation 
and subsequently the corresponding moment of inertia. 

Next, the normal modes are determined by an expansion of 
auxiliary characteristic functions (using N ev numbers! out of the N number 

Tit is probably more accurate not to use the last few auxiliary eigenfunctions 
which contain larger errors. 
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of 4*1 s ). A matrix eigenvalue problem results as given in Appendix II. The 

w 

computer program thus calculates various integrals entering the matrix 
equation. An eigenvalue subroutine yields the coefficients in series expan- 
sion of the normal modes (eigenvectors) and the corresponding frequency 
(eig envalue). 

Once the coefficients in the series expansion of the normal 
modes are determined, integrals containing these functions (Appendix II) can 
be evaluated to determine the parameters of the spring-mass mechanical 
model. This completes the main steps of the present program. 
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II. EXAMPLES 

Results of numerical examples are given in Appendix III. 

It is noted that for a spherical tank with a "folding 11 interface, the 
result is sensitive to the inputs, especially the interface points; piecewise 
uniform spacings of the interface points were used. With a 12 X 12 mesh, 
it takes about two minutes CDC 6600 time in the main program, TRIPOT. 
With 23 X 34 mesh, it would take about twenty minutes. 
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III. CONCLUSIONS 

The program seems to yield reasonably good results in the examples 
calculated with as few as 12 interface points. However, in some cases, it 
is sensitive to inputs which should be reasonably accurate. Applications to 
other cases remain to be examined. 
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IV. INSTRUCTIONS FOR USING COMPUTER PROGRAMS 
"TRIPOT 1 1 AND "SSHAPE" 

A. General 

TRIPOT is a program that computes the sloshing parameters for an 
axisymmetric tank. A finite -difference technique using a triangular mesh 
is used to solve the potential flow equations. The main restriction on the 
program is that the Bond number and the volume of contained liquid are both 
great enough so that the liquid covers the bottom of the tank. 

SSHAPE is a program used to compute the equilibrium free surface 
shape of the liquid. The output of SSHAPE is used as input for TRIPOT. 

B. TRIPOT 

A triangular mesh is laid over the liquid by a method to be described. 
Figures 4 a and 4 b, in Appendix III, show typical meshes for a cylindrical 
and a spheroidal tank. 

Both the liquid interface and the tank shape are given as discrete x, y 

t 

coordinates, x = 0, y = 0 being the point where the free surface intersects the tank 
centerline. The interface, arc length sj, is broken up into N points. The center - 
line, arclengths^, is broken up into M points. The bottom, S3, and side wall, s^, 
are specified at N and M points, respectively. For a spheroid, for example, in 
which there is no "bottom" or "side," the wall is arbitrarily divided into arc 
lengths S3ands4 such that s^/s^is approximately equal to s^/sg. 

For a "folded over" interface, additional accuracy can be obtained by 
using a non-uniform spacing of the N and M points. A coarse net is used 
on the interface, near the centerline and on the part of the wall not close 
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to the contact point, s^. Then, S£ and are divided similarly to Sj and s^, 
respectively. For example, 16 intervals (N = 17) can be used in SSHAPE to 
calculate the free surface shape. Then, the s, y coordinates of the 1, 3, 5, 
7, 9, 11, 12, 13, 14, 15, 16, 17 points can be used as input to TRIPOT with 
N = 12 (eleven intervals, 5 wide ones near the centerline and 6 narrow ones 
near the contact point). 

The interior grid points are generated by TRIPOT. In general, it is 
better to use near-equilateral triangles and avoid obtuse ones. Some trial- 
and-error work may be necessary to insure this. 


Subprograms Used 


SOREL - 


Solution of a system of equations by successive 
over -r elaxation 


SFIT 


Curve fitting by a second degree polynomial 


EIGEN - 


Eigenvalues and vectors of a real symmetric or non- 
symmetric matrix 


MATINV - Matrix inversion routine 

Tapes Used 

TAPE 6 - Restart write tape 
TAPE 7 - Restart read tape 



Less than 32, 000 case storage locations 



13 


Input Formats Used 


Card No. 


1 

2 

3t 

4t 

5t 

6f 

7 

8 
9 


F or mat 


(12A6) 

(415, 5F10.0) 
(8FI0. 0) 
(8F10. 0) 
(8F10. 0) 
(8F10. 0) 
(8F10. 0) 
(5E15. 8) 

(415, 5FI0. 0) 


C. TRIFOT- -Program Restart Procedure 

At seven different places in the program, results calculated thus far 
are saved on tape so that it is not necessary to calculate from the beginning 
in the event of catastrophic failure or insufficient alloted time. When using 
the restart feature it is necessary to interchange tape six with tape seven so 
that the tape which was written will become the tape being read. One method 
of effecting the interchange is by redefining LUN and KUN, which are the first 
two executable statements in the main program. During execution, the printed 
output will state the highest restart number which can be used. The input data 
card numbers to be used with each restart are given below. 

Restart Number Input Data Card Numbers 


1 

2 

3 

4 

5 

6 
7 


1,2, 7, 8, 9 
1 , 2 , 7 , 8, 9 

1 . 2 . 7 . 8. 9 
1 , 2 , 8 , 9 
1 , 2 , 8 , 9 

1 . 2 . 8. 9 
1 , 2 , 8, 9 


TSet of cards. 
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TABLE I. INPUT OF TRIPOT 

TRIPOT — Input Data Description 
Card FORTRAN Variable 


No. 

Symb ol 

Name 

Units 

Definition 

1 

ITITLE 



72 columns of identification 

2 

M 



Number of centerline points 


N 



Number of interface points 


ITA 



Maximum number of iterations 


IRS 



Restart number (initially zero) 


WA 


non- dim. 

Relaxation factor for linear 
a ppr oximati on s 


WB 


non -dim. 

Initial relaxation factor for 
nonlinear solution 


EPS 


non- dim. 

Convergence factor 


RLGTH 


length 

Reference length used for non- 
dimensionalization, e. g. , 
maximum tank radius 


BOND 

n b 

non -dim. 

Bond number 

3t 

X(1,J), Y(1,J) 

( x , y)j=i» n 

length 

Coordinates for interface^; 
j = 1 at tank center, increasing 
outward 

4t 

X(I,N), Y(I, N) 

(x, y)i=2, M-l 

length 

Coordinates for wall; i = 1 at 
interface increasing downward 

5t 

X(M, J), Y(M,J) 

( x > y)i =1 *N 

length 

Coordinates for centerline, j = 1 
at tank center, increasing outward 

6t 

X(I,1), Y(I, 1) 

(x, y)i=2, M-l 

length 

Coordinates for centerline, i = 1 


at interface increasing downward 


tSet of cards. 

;|It is recommended to use 5 digits (or more) obtained from output of interface 
program. Wall coordinates are not as critical. 
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TABLE I (cont'd) 


Card 

No. 

FORTRAN 

Symbol 

Variable 

Name 

Units 

Definition 

7 

GFLC 


non -dim. 

= 0, compute T ; input dummy F 
= 1, input T from SSHAPE out- 
put, or if interface is flat input 

r = o. 


GAMMA 

r 

non-dim. 

Contact point constant 


THETAC 

9c 

degrees 

Contact angle; must be greater 
than Z° 

8 

ZCG 


length 

Center of gravity of liquid f 


YOL 


length 

cube 

Volume of liquid f 


DEN 


mass /length 
cube 

Density of liquid 


RHOU 


mass /length 
cube 

Density of ullage fluid 

9 

t 

NEV 


non -dim. 

Number of eigenvectors used 


f Obtained as output from SSHAPE. 
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D. TRIPOT Input Example and General Instructions 

Example fox a Spheroidal Tank and General Inputf Instructions are 
given below: 


1. TEST CASE xx, x Full, Bond number = x (see input description) 


2. 

a. 

M = 9 

no. of side wall points, see 5 below 


b. 

N = 12 

no. of interface points selected, see 4 below; 
usually greater than 10 


c. 

IT A = 200 

suggested value; usually sufficient 


d. 

IRS = 0 

always zero unless using restart procedure 


e. 

WA = 1.8 

suggested value 


f. 

WB = 1.3 

suggested value 


g- 

EPS = 1 x 10“ 6 

suggested value 


h. 

RLGTH =1.0 

(see input data description) 


i. 

Bond = 5 

specified value 

3. 

X(1,J), Y(l, J) 

selected interface points from output of SSHAPE 

4. 

X(I,N), Y(I,N) 

side wall points selected after interface is 
calculated; a plot of interface, tank wall and 


center is desirable for selecting wall points. 
Suggest division near contact point be nearly 
an isocicle or form an acute triangular. It 
is not necessary to have precisely equal arc 
divisions. Usually, the meshes may be finer 
'near the contact point. The points can be read 
off from curve (or calculated from analytic 
expression); accuracy is not as critical as the 
interface points. 


t TRIPOT input data descriptions are given in Table I. 
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5. X(M, J), Y(M, J) bottom wall points; the ratio of bottom wall to 

side wall m arc length is approximately that of 
interface to wetted centerline (eye observation 
is usually sufficient). The divisions on the 
bottom are suggested to be approximately 
parallel to that on the interface. Points can 
usually be read off from curve plotted. 


6. 

X(I, 

1), Y(I , 1 ) 

wetted centerline points. The divisions on 
the centerline is suggested to be approximately 
parallel to that on the side wall. 

7. 

a. 

GFLC = 1.0 

(see input data description) unity using output 
from SSHAPE 


b. 

GAMMA = -16. 12 

output CGAM from SSHAPE 


c. 

THETAC = 5 ° 

specified contact angle (no less than 2°) 

8. 

a. 

ZCG = -. 01 

output from SSHAPE; c. g. of liquid, negative 
below center of interface 


b. 

VOL = 1.35 

output from SSHAPE; liquid volume 


c. 

DEN = 1.0 

if RHOU = 0, the density of the liquid DEN 
can be set to 1 without affecting result 


d. 

RHOU = 0 

density of ullage fluid usually neglected (zero) 

9. 

NEY 

= 7 

number of eigenvectors used; for N > 11, 
N ev = 7 has been selected. For partial 


convergence test on a finer mesh N ev = 7 
can be retained and for the convergence 
test on Galerkin procedure, N ev may be 
increased. However, it may be advis- 
able to use N ev < N-3 or smaller to avoid 
large inaccuracy m the higher auxiliary 
characteristic functions. 


E. TRIPOT Output 

TRIPOT --Program Printed Output 
1. Input Data - TITTLE, M, N, ITA, WA, WB, EPS, IRS 


2 . 


(x, y) coordinates for linear approximations to the mesh 



3. 
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(x, y) coordinates for nonlinear solution to the mesh 

4. Input Data (compute option) - GAMMA 

5. Eigenvalues X. and eigenvectors 4^. on V 

•* , j 

6. Input Data - BOND, ZCG, VOL, DEN, RHOU 

2 

7. Eigenvalues' S2 2 = and e ig envect0 rs C, 

g 3 

2 

8. SIF, I F /M F h o 

9. SIO, I 0 /M F h Q 2 

10. MK/MFj - nondimensional i*k slosh mass, M^/M^ 

11. ZK , slosh mass location, z with i = k = 1, 2. . . NEV 

1 K 

12. Z0, location of rigid mass, 
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APPENDIX I 

A COMPUTER PROGRAM FOR GENERATING THE INTERFACE 
SHAPE AND ASSOCIATED QUANTITIES 
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A. Governing (Mean) Interface Equation 

The interface differential equation is given byf 


,,„ 2y2 + 3(y») 2 y> 


cot (0) [ y2 + {y 1 ) ] 


+ I 

y 


Bn (y cos (0) - y Q ) - 


2k 


0 


ly* + (y» ) 2 ] 3/2 (1-1) 


Converting from second order to first order, 

let 


Yi =y 

y 2 = y 


I 


then 


y[ - y' 

y'z = y" 


(1-2) 


or 




2y 2 + 3y| 

yi 


cot 

n 


(0) ly\ + y| 1 


(1-3) 


+ 





( 6 ) - y 0 ^ 



( y 2 + y |) 3 / 2 


tHastings, L. J. and Rutherford, R. , III, "Low Gravity Liquid- Vapor Inter- 
face Shapes in Axisymmetric Containers and a Computer Solution, 11 NASA 
Technical Memorandum, NASA TM X-53790, October 7, 1968. In the equa- 
tion to follow, y is the radial distance from the center of the top of the tank, 

0 is the angle counterclockwise from the tank centerline with origin at the 
center of the top. 
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The solution of the first order differential equations (1-2) and (1-3) is 
obtained by using the Runge-Kutta-Gill fourth order method; the boundary 
of the tank is approximated by a polygon. In general, a new polygonal 
approximation is required below the interface according to the description 
aforementioned in the main text. 

B. Program Notes 

In addition to the main program, the following subprograms were 

used: 

RKLDEQ - Runge-Kutta-Gill differential equation solver 
SFIT - curve fitting by a 2nd degree polynomial. 

The input data description of program SSHAPE is given in Table II and the 
output of program SSHAPE is given in Table III, respectively. Some of the 
program symbols which were not given previously are defined in Table IV 
and an example keypunch form is given in Table V. 

C. Input Instructions for SSHAPE 

The input for a spherical tank and general instructions are given 

below: 


1. 

NP = 37 

Total number of points specifying the tank 
starting from center of bottom to center of 
top; varies with each case 

2. 

YBX. 

Horizontal and vertical coordinates of tank; 


1 

should be a close polygonal approximation to 

3. 

YBZ, 

the tank, 1 = 1 starting from center of tank 
bottom 

4. 

a. a = 5° 

Contact angle specified 


b. KO = 0 

Usually zero 
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c. 

YO = 1 . 0 

Initial guess at distance between center of 
interface and center of tank top 

d. 

BN = 5 

Bond number specified 

e. 

BETAD = 5/8 = . 625 

Empty fraction specified 

f. 

DKO = 0.2 

Suggested 

g- 

DBC = . 005 

Convergent criterion 

h. 

DYO = . 05 

Suggested 

5. a. 

THETA = 0 

Usually zero 

b. 

DTHETA = . 05 

Suggested 

c. 

RLGTH = 1 

(See input data description) 

6. a. 

N = 2 

Always 

b. 

IBOPT = 1 

Unity if calculating empty fraction (see 
input data description) 

7. a. 

NN{1) = 17 

(Only 12 of these 17 were selected for 
input in TRIPOT) 

b. 

NN( 2 ) = 33 

For a finer mesh if desired 

D. 

Listing and Sample Input 



These are included, after those of TRIPOT at the end of the report. 
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TABLE II 

INPUT DATA DESCRIPTION OF PROGRAM SSHAPE 


Card 

No. 

FORTRAN 

Symbol 

Variable 

Name 

Units 

Definition 

1 

NP 



No. of points on container boundary 

2 

YBX 

l 


length 

x-coordinate container boundary, 
positive outward from top center, 
i = 1 at bottom center 

3 

YBZ 1 


length 

z -coordinate container boundary, 
positive downward from top center, 
i = I at bottom center 

4 

ALPHA 

a 

degrees 

Desired contact angle, 9 C , must 
be greater than 2° 


KO 

K o 

non -dim. 

Assumed initial guess of parameter 
related to curvature k Q at inter- 
face center point, i. e. k Q y 0 . 
Usually zero can be used. 


YO 

Yo 

length 

Estimated distance from center of 
tank top to interface centerpoint 


BN 


non-dim. 

Bond number, N 

B e 


BETAD 

Pd 

non-dim. 

Desired empty fraction 


DKO 

> 

0 

non -dim. 

Increment for K Q , a typical value 
is 0. 2. 


DBC 

e pD 

non -dim. 

Convergence criterion for (3^. 
{A typical value is 0. 005) 


DYO 

*Yo 

non-dim. 

Increment for y Q . {A typical 

value is 0. 05 R ). 

o' 

5 

THETA 

0 

degrees 

Initial angle measured from vertical 
axis to y . Usually zero. 
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TABLE II 

INPUT DATA DESCRIPTION OF PROGRAM SSHAPE (Cont’d) 


Card 

No. 

FORTRAN 

Symbol 

Variable 

Name 

Units 

Definition 


DTHETA 

A0 

degrees 

Increment for 0. (A typical 
value is . 05) 


RLGTH 

R o 

length 

Characteristic container dimension 

6 

N 

n 


No. of equations 


IBOPT 



Option for calculation of empty 
fraction P . If = 0, suppress 
calculation"!"; if = 0, calculate P 

7 

NN(1) 

NN(2) 



Two choices of no. of interface 
points of equal intervals of which 
the locations are desired. 


tin case y 0 is known from either another theory or experiments, one may proceed 
the calculation with a dummy mass fraction. In theoretical predictions, usually desired 
empty fraction is given, then we must use IBOPT = 1 to calculate p for each assumed 
y Q in the program, while the input y Q is an approximate guess. 
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TABLE III. OUTPUT OF PROGRAM SSHAPE 


Printed Output 


1 . 

Input data ALPHA, KO, YO, BN, THETA, 

Y(l), 

Y(2) 

2. 

THETA, XOR, ZOR, Y{1), Y(2) 




3. 

GAMMA, PHI, ALPHA, YO, KO 




4. 

Extrapolated YO (non-dimensional) 




5. 

VU1, VU2, VU, VL, VT 




6. 

BETA 




7. 

G, RS, S, XR, XY 




8. 

K2F, FRR2 (interface curvature n 
point) 

and 

d 2 F 

d R 2 

on F at the contact 

9 

K2B, ZRR2 (wall curvature and 

CGAM ( r ) 

d 2 Z 

on 

wall at the contact point) 

10. 

V 


11. 

Dimensional YO 




12. 

XRB, XYB 




13. 

ZCGT, BZCGU, ZCG 





FORTRAN Variable 

Symbol Name Units Description 


S 

s 

length 

l 



G. 

G 

length 

1 


square 

RS 

i 

R s 

length 


Arc length, i = 1 at interface 
center point 

Sum of curvatures squared, 1 = 1 
at interface center point 

dR/dS, i = 1 at interface center 
point 
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TABLE in. OUTPUT OF PROGRAM SSHAPE (Cont’d) 


FORTRAN 

Symbol 

Variable 

Name 

Units 

Description 

CGAM 

r 

non -dim. 

Contact point constant 

XR 

XR 

length 

Horizontal distance from container 
vertical axis to liquid -vapor inter- 
face with origin at interface center 
point 

XY 

XY 

length 

Vertical distance from x-axis to 
liquid vapor interface with origin 
at interface center point 

XRB 

XRB 

length 

x-coordinates of container boundary 
with origin at interface center point 

XYB 

XYB 

length 

y-coordmate of container boundary 
with origin at interface center point 

ZCGT 

Z 

CGT 

length 

Container c. g. 

BZCGU 

P Z CGU 

length 

Mass fraction times distance of 
ullage c.g. upward from interface 
center point 

ZCG 

Z 

CG 

length 

Distance of liquid c. g. upward 
positive from interface center 
point 
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TABLE IV. DEFINITION OF TERMS IN SSHAPE 
Some of the program symbols which were not defined previously. 


FORTRAN 

Symbol 

Variable 

Name 

' Units 

Definition ' 

XOR 

X/R 

0 

non-dim. 

Horizontal distance from con- 
tainer vertical axis to interface 
point 

ZOR 

Z/R 

o 

non -dim. 

Vertical distance (downward 
positive) from x-axis to interface 
point, origin at top center. 

Y(l) 

y 

non -dim. 

Radial distance from top center 
to interface point 

Y(2) 

¥' 

non -dim. 

First derivative of y with 
respect to 9 

GAMMA 

y 

degrees 

Angle measured from vertical 
and tangent to y(0) at contact 

point 

PHI 

4> 

degrees 

Angle measured from vertical 
and tangent to y (0) at con- 
tact point ® 

VU1 


length 

cube 

Vapor volume to point of contact 

VU2 


length 

cube 

Vapor volume from point of con- 
tact to boundary 

vu 


length 

cube 

Vapor volume 

VL 


length * 
cube 

Liquid volume 

VT 


length 

cube 

Total container volume 



TABLE V. EXAMPLE KEYPUNCH FORM FOR SSHAPE 


COMPUTATIONS LABORATORY 


1 ■■ 

PROBLEM Ex-siriplp Keypunch Fo-rvn 

PAGE OF 

PROGRAMMER 

DATE Q f 2 0 / 



I 

I 


O 
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APPENDIX II 

A THEORYt FOR LOW-GRAVITY FUEL SLOSHING IN AN 
ARBITRARY AXISYMMETRIC RIGID TANK 

Wen-Hwa Chu 


fThis theory and mam results have already been published as ASME 
paper 70-APM-EEE to appear in the Journal of Applied Mechanics and con- 
tains major extensions and modifications of Technical Report No. 8. It is 
included for convenience of the sponsor and the readers of this report. 
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Introduction 

Thk behavior and consequences of fuel sloshing in 
rochets under a high effective gravity were recognised problems 
which have been quite well understood [ 1-81 1 The pioblem of 
low-giavifcy fuel sloshing, characterized by the significant role of 
interfacial tension, is now a subject of importance for application 
to coasting rockets or orbital stations 
The equilibrium behavior of fluids at zero and/or low giavity 
has been studied m references 14-7] The theoietical determina- 
tion of an equilibrium interface shape is nonhneai and lequires a 
trial and erioi procedure for a given contact angle [">, 6] 
Satterlee and Reynolds [8] have successfully solved the free 
sloshing pioblem in cylindncal containers undei low gravity and 
foimulated a variational principle for this puipose Yeh [9], 
using a similar approach, solved the fiee and foited sloshing 
problem undei low-gravity conditions, without force and moment 
or an equivalent mechanical model Dodge and Garza [10, 11] 
pet formed force measurements undei simulated low-giavity con- 
ditions and predicted forces and moment for circular cyhnducal 
tanks undei lateral (translational) motion The equivalent 
spring-mass model was given in [10] Additional work by 
Dodge and Gaiza for other special tanks was given in [12, lol A 
finite-difference approach with application to a henusphencally 
bottomed cylindrical tank and spheroidal tanks was given by 
Concus, Crane, and Satterlee in [14, 15] 

These investigations indicate a need of a progi am for a general 
a\isymmetric tank A preliminary study on liquid sloshing in 
an arbitrary axisymmetuc tank was reported in [16], but it is 


Governing Equations 

Assuming inotatunial incompressible flow, theie is a space- 
fixed velocity potential 1 ip satisfying the Laplace equation in bot h 
space-fixed and tank-fixed cool dmates 

V,-<£ = 0, V’</> = 0, 

_ a* a 2 a* „ a* a* a 5 

(7l — 1 J V 1 1 + — - ( I 1 

* bx,' by.' dz, 2 ’ as* by- bz 1 

As m thin airfoil theory, the velocity potential can be obtained 
by imposing boundary conditions on the initial or mean position, 
but the hydrostatu pressme due to gravity possesses components 
along both the tank axis z and the lateial axis x, Fig 1, for pitch- 
ing oscillations The hneai i/ed Bernoulli's equation states 

bi> ~ 

p - pi + p ~ + pg(z - x6„) = 0 <2) 5 


and 


P - P«i + P» + p„Sf(z - xd„) = 0 


O) 


for the liquid and the ullage, respectively, and pi, p„| aie c (in- 
stants 


Boundary Conditions 

The linearized mteifaoe kinematic condition states that 
bh ^ b<j> 
bt bn 


on F 


<4)‘ 



The mtei fai e dj namic condition states that 

— p~ + p+ — (JK = OKo + ok' on F, (5) 

1 

For the “mean” interface location/ (in general, pi = pi° + pi', 
pi°, pi' being constants), 

| — <rxo -f- (p — pjff/ — (pi° - p„i°) = 0 on F 16) 

wheie the ( in vat me of the mean interface, x 0 , is axisymmetric 
and 


1 b 
x'o = - — < - 


b> 




r &$ \5s as 2 a$ a* 2 / 


Equation (G) holds for r — 0, thus 

The hneai t/ed interface dynamic condition is then 

dtp dtbu di 

-(?«' - P-l') - + P T7 - pu — + (p - P«)gh — 

Ol ol O S 


limited to 1 1 Auslat lonal oscillation^ It is the object of the 
pieseiU papei to piesent a semi numerical approach for an arbi- 
tral y axisymmetuc tank with simplified foice and moment cal- 
culations and the i exultant mechanical model foi both pitching 
and tianslational oscillations A general computer program 
utilizing Winslow method (17] will be completed to obtain slosh- 
ing frequencies, slosh mass, and mass-height, foi which a brief 
description is given in Appendix A 

1 Numbers in brackets designate References at end of paper 

Contributed by the Applied Mechanics Division for publication 

(without presentation) in the Journal of Applied Mechanics 


— (p - p*)gx6 t = 0 (7)* 


3 For example, gives xeloutj toinponent in x-direction with 

respect to N]> t ue-fi\ed coordm ites 

* Chu, W H , “Free Surface Condition for Sloshing Resulting 
From Pitching and Some Corrections,” ARS Journal , Vol 30, Nov 
i960, p]> 1093-1094 

4 A vertical interface description was used but is onh successful 
when (dF/dR)~* does not vanish on the interface 

5 For sinusoidal oscillations and m — 1, h = 0, = <£« — 0, x — 0, 

and K f = 0 at point I, thus pi' — j> WI ' = 0 For other m values, 
— pt* 4- puj° = tTKi\ which will be omitted until needed 
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wheie the peitiubation curvatme foi cos (mO) vauatiou is to the 
first order 


5 2 X 1 5r 5X 1 _ OTi 

5s 2 r 5s 5s ^ r 2 50 2 


+ 


[(;!)’ + (I 


faf 

5s 2 


VWrVl 

5s 5s 1 / J 


h 


(7o) £ 


m being unity foi lateial excitation of a ngid tank At. point I, 
the ongm, h = 0, <t> = 0, x' = 0, and thus pi = p UI Foi most 
analyses, p„ — 0 was assumed We shall assume the impulsive 
pressure in the tillage is negligible, i e , = 0 Then, foi sinu- 

soidal oscillations with h, <f> piopmtional to cos (cot) and sm (cot), 
lespectively, equations (7), (7a), and (4) yield 


(I 

R bs \ 3s / 


-w> ff+GH 


4- Nb,H 


, 5 R 

5s 


-j- 8 y N a Jl cos 8 + = 0 

where 


(7&) 2 


/ 1 5F\ 2 /5K &F _ 5F 5^y 

V? 5s / + \5s 5s 5 5s 5s 2 / 


jl , r y 

It'- (1 + F„ 2 ) + L(1 + W /2 J 


(7c) 


It is noted that, foi lateral oscillations, tn — 1, solutions are pio- 
portion al to cos 6 

The boundary condition on the wall is that the relative normal 

3x , , , 3z 

velocity be /eio, 1 e , with cos (n, x) — ^ and cos \n,z) — ^ 


and 


b(j> 3x 
3n ° 


(S) 


3c/> 

3n 



3x 

bn 


bz y 
bn) 


(9) 


for translational and pitching oscillations, respectively 

In addition, there is an interface contact point condition which 
takes the form (8, 9, and 15] 

6 The instantaneous interface can be described by the position vec- 
tor r/a [R - HF<\ cos 9i+ [R- HF S ] sm 0} + [F -I- R* - 

dR/ds, F, « dR/ds The curvatures can then be derived from the 
fundamental magnitudes (18] neglecting higher-order terms The 
linearized total curvature is given hy this equation, uhich wts first 
derived in [151 by a_chfferent procedure 

» For sinusoidal oscillations, without loss of generality, zo, 0 V , 
are assumed to be proportional to sin (orf) while h is proportion il to 

COS («0 


— Nomenclature 

a = teference length, say, maxi- 
mum radius of tank 
J 1 = rdrdO 
if A = d.i /a- (a -scalai ) 
tIS = i-D surface element, eg, 
rddds 

(IS — dS/a 2 , nondimensional siw- 
face element (a scalar) 

F = equihbuum (mean) mteifaie 
of f/a 

F t — instantaneous mteifaie 
F„ = hoi i/on t a! foice defined by 
equation (19) 

h „ = UlF/dlt), slope of F in tbe 
genei atnx plane 

t x = x-componenl of foice on tank 
f = equihluium (mean) mteifaie 
elevation measuied along 
vei tical axis 

g = giavitationai aiceleiation 
H — amplitude of h/a, nondimen- 
sional slosh height 
h = mteifaie peituibation noimal 
to equihbi mm interface 
h a = lefeienie length, say', depth 
of liquid at eentei of tank 
il la = i igid ma-ss and moment of in- 
ertia of mechanic al model 
1 I f = liquid mass 
M x = pitching moment about i/-axis 
m t = It h slosh mass 
Xa = Bond mimbei, ptFg/cx 
n = outei noi mal 
Xu = n/a, nondimeiisioital noimal 
distant e 
p = pressuie 

pt = equilibrium liquid piessiue at 
origin — a constant 
p u = ullage prcssme 


= equihbuum ullage piessiue at 
origin — a constant 
It - i /a, nondimensional lidms 
i, 0, z = tank fixed cylindrical coordi- 
nates 

i = arc length nondimensionah/cd 
by a 

s = ate length (a scalar) 
f = time 

V = volume of liquid divided by a 2 
V r = liquid volume 
ir = wall wetted by liquid 
IF, - instantaneous wetted wall be- 
low instantaneous intei face, 
F, 

to = tianslational amplitude in x,- 
dncction 

„g„z, - space-fixed reitaugulai eooi- 
dmates 

I’ = yu, nondimensional contact 
point constant 

y = hysteresis coefficient oi con- 
tact point constant 
A p = density diffeience, p — p„ 

S„ = Kroriecker delta 
6, = sign of a 2 , coo (n, z),oi Dz/dn 
0„ = amplitude of pitching about 

i;-a\is 

x = mean imvatuie 
x' = perturbation of mean imva- 
tme 

X, = yth eigenvalue (m = l) 

\ mj = jth eigenvalue couesponds to 
with cneumfeiential mode 
p = liquid density 
P„ = density of ullage Huid (vapm 
oi gas; 

er = siu face tension 
d 1 = amplitude of nondimensional 
velocity potential, <j>/wa 2 


= amplitude of nondimensional 
potential 4>il b,a> 

(f> iV - see equation (15) 

= amplitude of nondimensional 
potential <j>°/ua 2 
4> — velocity potential 
fa = velocity' potential of the 4lh 
natural mode 

fa = additional velocity potential 
due to interface movement 
fa = velocity potential of liquid 
with a frozen intei face 
xp — velocity potential of auxihaiy 
eigenfunctions 

xj/ / = jth auxiliary charade! istii 
function 

(J 2 = pa'dr/a, product of Bond 
number, pa 1 g/<r, and fie- 
quency' parametei, a Fa/g 
ip = fiequency of oscillation 
c oi_ = 4th natural fiequency 

Subser/pfs 

( ), = ( ) at veitex of equihbuum 

mteiface (ougm) 

( )n = ( ) at contact point in genei - 

atnx plane 

( )ci = ( ) related to center of grav- 

ity 

( ), = effective value of ( ) 

( )x = ( )o « F 

( ) m = ( ) associated with cos (mo) 

mode 

( ) ji = ( ) related to pitching 

( ) r = ( ) related to translation 

( )n = ( )onlt' 

( )« = ( ) related to ullage 

( )_ = ( ) just below interface 

( )+ = ( ) just above mteiface 
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where y may be a frequency-dependent constant Howevei, if 
the contact angle remains constant, w e can show that (see Appen- 
dix C 8 ) 


sm 9 C { 


a + z T -)' h 


— cos 0. 


(i + jV)‘ / ’ In] 


r = ya 


(10 a y 

(105) 


whei e 


dr’ h dr’ " 


dh 

dr 2 ’ 


;,c 


dr 2 ’ 


z on W, f on P 


Method of Solution 

We shall decompose tj> into two parts, <f> ' and <j>° <f>° is the 

velocity potential corresponding to a liquid contained by a rigid 
mean intei face and the tank walls Therefore, it satisfies the 
Laplace equation and the nonhomogeneous boundary condition 
on the contour, equation (8) for translation and equation (9) for 
pitching on F c and W e It is noted that 

(f> T ° — x<& (II) 

while <£ p ° can be constructed numerically 

<j>' is the perturbed velocity potential due to sloshing which is 
governed by the Laplace equation, the zero normal velocity con- 
dition at the wall and the resultant interface condition The 
fiist two are satisfied by expansion in normal modes and the last 
by equations (15) and (16) In this paper, the normal modes are 
determined by expansion into a set of auxiliary characteristic 
functions 10 i !/ v which is orthogonal on the curved interface, and 
satisfies the Laplace equation and the zero normal velocity wall 
condition 11 The intei face condition governing normal modes is 
second order [see equation (76)] subject to zero H at center and 
the contact point condition, equation (10a) at wall The former 
is satisfied as ^ vanishes along tank center line The interface 
equation is approximately satisfied by a modified Galerkin 
method (the strict Galerkin method is given m [19)), which is 
illustrated m Appendix B In this method, the contact point 
condition was imposed on the whole series, not term by term 
The numencal results shown later substantiate the method used 

The force and moment are obtained by integration of pressure, 
not only on the wall, but also on the intei face since the direct sur- 
face tension force and moment on the tank is equivalent to those 
on the interface due to pressure, assuming the interface inertia is 
negligible, as well as the interface mass To put results m the 
mechanical model form, the divergence theorem has been most 
useful (with some manipulations) 


Analytical Results 

Free Oscillations. For free oscillations, the natural mode 0* 
expanded into a truncated senes of the auxiliary eigenfunctions, 
le, 


8 Appendix is not in order of mention 

8 An equivalent form w T as first derived in reference [15 J 

10 For direct application of the Winslow method [17J, we impose the 
simpler normal derivation condition, di/'j/bno = on F and used 
the well-known influence coefficient technique to determine the eigen- 
vector on the mtersurface, the eigenvalue A, 

11 Strictly speaking, the boundary conditions are imposed in Wins- 
low method, white the solution not proven to satisfy these conditions 
is correct on physical grounds (see Appendix A) 


h Jmx 

Hh = — = - ^ cos 

° 3=1 

ci, is the fcth eigenvector of the following matrix equation obtained 
by a modified Galerkin method (Appendix B) from integrating 
the nondomensional equation (76) with 8„ = 0 and weighting 
function if/njlS/am, 2 


+ t7„„] + m*[e w J 

d ~ Aplfrnul ^*l^m tf l | (Cj} — 0, 1 , J — 1 to Jmx (13) 


where 







'f'mjl'mit l 

VTTf? 


dS, €i 





V mt 


27rX m „ 

Qm, 8 


[R'l'mrfm,] 


(13a) 1! 

(136) 

(13c) 


Ym,, 


Ci m , 


•if. 




J (13d) 



(13e) 



t/'m, 8 dS 


and m — 1 for lateral excitation of a rigid tank 


(13/) 


Forced Oscillations Let 

K m x 3 os J mx 

V = -“ a2 E d &> I 7w = E CL d> 

A = 1 “ 3 = 1 

(14a, b, c ) 

ra order to satisfy the interface condition that 

An Y 

E - C 8 )#* = + ftiV a — - 9 V = (15) 

k=l p a 


e» = 0 for translational oscillation, e, = 1 foi pitching oscillation 
We have by the Galerkin procedute [19, equation (115 5), p 435] 


Kmz /* 

E 4 

t. - 1 j f 


= f ‘Pxi-HJdS, l = 1 to K mz (16) 
F J F 


A weighting function of Hi was used taking advantage of the 
biorthogonal relation, equation (17) This assures the first K„ z 
components of the error as senes in <pi be zero d k can he solved 
from equation (16) by matrix inversion There is no need of 
storing information of i/q inside the fluid domain as only the force 
and moment are of interest It is noted in the limit [8 and 9j 


L 


$,(-ff,)dS = s lt 


J>- 


H k )dS 


(17> 13 


11 The orthogonality property of thus, can be easily proved 
[16) as in the high-(7 ease cos 1 (7rt0) in the integrand are uniformly 
omitted as they contribute the same factor, J /i 
18 This is referred to as biorthogonal relation 
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then, 


d,= 


L 

L 


ZA-HJdS 


^(-H t )dS 


(18) 


which w at. utilized in proving that a unique spring-mass system 
exists for both pitching and tianslation 
Force and Moment The force and moment exerted by a spring- 
mass system, Fig 2, w ithout damping were w ntten in the follow- 
ing foim [20] and remain valid if contiibutions due to duect sur- 
face tension me included 


do I F ( fflo Zo . nifc z k \ (26) 

MM ~ M P hJ ~ \M k ho- A Vi M e ho'-) 


The force due to liquid pressure and direct surface tension u> 


K 


J, 


W'+Ft 


Sa: 


(27)“ 


and the moment due to liquid pressure and direct surface tension 
is 



(28) 



it can be shown that 


where 


3. _ p I 
UrfL V 


(29) 


(29o) 


3tr = ~ f $ t ~d$3Z ~ J f c*,A f -^dS (296)“ 
PS J ir+F P*%ti Jf “ 

d, PS = j $*(- 


v = V,J(,a\ 


and 


where 


1* = J_ (iA 

a m k \Y ) 


-H k )d A (29c, d) 


(30) 


M p 


k’ = d k t X! 

j = l 

J , ( z bx x bz\ 

iA~T r- ) dS 

F+W \ a dn a hn J 


(30a) 

(306) 


and that 


1 C - / z br x bz \ 

I F = M F a'-I t \ I F * S ? j V ^ 


(31a, 6) 


In deriving the mechanical model, p„ has been set to zero A 
simple modification can be made for small ullage density by using 

Ap 

Nb,= —X b 

P 

the effective Bond number instead of the Bond number based on 
the density of the liquid, provided that the dynamic pressure due 


with rigid mass mo, its location z 0 , and moment of 

by 

mp . zn k 

Mr ~ Al Mf 

zo^iTzco,. ft m k ~ 1 

ho ?Ho L ('0 ho Mf J 


inertia Id given 


(24) 


(25) 


M f 


11 It is important to note that great simplification in results, when 
reduced to the mechanical model for both translation and pitching, is 
achieved by consideration of balance of forces acting on the thin inter- 
face as a free body The force and moment on the tank due to direct 
surface tension must be eciual to those acting on the interface by 
liquid pressure These were not included m reference [15) 

15 For finite Jm„ it was found that dir, determined by matrix in- 
version of equation (16) without using biorthogonal relation, vields 
results in better agreement with Dodge’s theory [12) than equation 
(296) which is correct in the limit 
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to ullage motion is negligible 

It is noted that to prove the mechanical model for pitching 
motion fiequent application of divergence theorem combined with 
diffeiential pioperties of the auxiliary characteristic functions and 
the coordinate functions x and z were used For brevity, the 
details are not included in this paper 

Numerical Examples The present computer program has been 
checked out for several examples m Appendix III. These include 
cylindrical tanks, spherical tanks, and spheroidal tanks Good 
agreement of first spring-mass with known theory or experiments 
was obtained with as few as 11 or 12 interface net points and the 
first seven auxiliary characteristic functions Finer nets may 
>ield higher accuracies, especially for the higher modes 

Machine Tima For a 12 X 18 mesh and a 12 X 12 mesh, the 
CDC-6600 central process time is a little over 2 mm, while for 
the 23 X 34 mesh it is about 21 nun Most of the computing 
time was expended for the generation of influence coefficients, 
each of which is a Neumann problem However, the influence 
coefficient method may be more convenient than the inversion of 
a laige matrix, if not faster No computer running time was 
reported in references [14 and lo], which finds natural modes by 
finite-differences arid (partial) matrix inversion 

Conclusion 

It seems that the pie^ent method yielded a practical way of 
computing the fundamental natuial fiequency, the first slosh 
mass, and its location Highei masses and locations are usually 
not needed foi design purpose* and can be obtained by using finer 
meshes and longer machine time A computei program utilizing- 
tnangulai meshes and Winslow method [171 has been success- 
fully employed and is expected to be completed in the near future 
foi the titled problem Howevei, the piesent logical diagram 
may be limited to a convex a\i:>yminetiic tank for good accuracy 
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APPENDIX A 

Brief Description of a Computer Program 

The following steps of a computei piogiam aie bnefh de- 
scribed 

Construction of a Triangular Mesh The tnangulai mesh K genet - 
ated as described m leference [17] eveept that a simple pmallelo- 
gram is used as the logical diagiam, Fig 3 Foi a c:\lnidi n.al 
tank at Bond number 100, the physical diagiam is shown in Fig 


LOGICAL PLAICE 



Fig 3 Simple logical diagram for triangular mesh 
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4(a) For a sphetoidal tank (c = 0 3) at Bond number .'>, a ln- 
angular mesh is shown in Fig 4(6) The lengths of the edge of 
the paiallelogram can be adjusted for each individual case to yield 
desned tnangular meshes A continuous wall need-, to be 
broken into two parts for the logical diagiam This only affects 
the local distubution of the tuangulai mesh and has shown to 
yield good results fot a spheioidal tank, as well as a cylnidncal 
tank 

Construction of Auxiliary Characteristic Functions The chaiacteils- 
ttc functions <p satisfy 


angle cannot be constructed graphically but results of decrease 
mesh sue give closer and closer approximations to the intei face 
and would probably lead to the correct limiting value 

For an interior joint, tj, If/ = if/,,, f/ L = Wh }), n = r tb, ]), 
r = r.jl 

6 

2 o) t (ft - f>) = 0 (35) 

L=l T 'J 

where 


<1 

II 

o 

(32) 

r 1 - = 0 on W 
Olio 

(33) 

- — = Ai p on F 
orio 

(34) 


f/ can be solved numerically with the constructed triangular 
mesh by Winslow method [17] Contact point is treated as one 
of the mesh points, as are the other boundary points 1S Hence, 
bf//bn may be discontinuous at the contact point Zero contact 



-2 4J 

Fig 4(a) Physical diagram of triangular mesh— cylindrical tank 


A tJ = area of the tjih dodecagon [17] inside the fluid domain 
r tJ = radius of the tjth. point 

wjl = £ (At,?*, cot 6k + \k-\fk-\ cot <Tjl) A- — 1 to G (35a) 17 

ft - i (r„ + n 4- r t+ i) Xi = 1 (356, c) 

0t, ert, Fig 3, can be expressed in terms of li, wt+i, Il-i, m-i, 
and st 

For interface point, 

£ «*(& A < $ 
r = i r '> 

+ (S),., [ 5 * + 5' <“> 

where 

Xe = Xi — X* = 0 , Xa ~ X 4 — Xs = 1 

To solve for the eigenfunctions on the interface, we use the in- 
fluence coefficient method in which {bf//bn) 1 % = 0 except {inf// 
bn)u, = X for the yth column of the influence matrix A stan- 
dard eigenvalue problem involving only the interface points, ex- 
cluding xpi 1 at r = 0, is needed to obtain the eigenvalues X, and 
eigenvectors f 1 , Knowing the jth eigenvector on the interface, 
the conesponding value of i/i, on the wall can be easily solved 
numerically again by the method of overrelaxntion 
For ij tli point on the tank wall 

ft m2 

2 utf.'l't - i>) A,,'!' - 0 (37) 

i=i T <i 


A3 = — A„ = 0 and Ai = A2 = A e = 1 on the bottom wall 

Xi = As = A6 — 0 and Ai = A? = X 3 = 1 on the side wall 


On center line, r — 0, 


^ = 0 for m > 1 
dxf/ 

— = 0 for 7/i=0 
or 


Ai — A* — A* — 0 and A* = A,> — Ag = 1 


(38) 


At contact point t = 1, j = 

E - *> - ^ A v * + Q -) r„ = 0 (39) 

Xs — 1, Xi = X- = Xi — X* — Xe — 0 


Calculation of Interface Shape A program based on the theory of 
reference [6] was written to generate the interface shape with de- 
sired net points for a given empty fraction or central depth and a 

16 This represents conservation of mass in the tn ingle cont lining 
contact point Experience shows that the t ont ict point condition 
should not be imposed on <//, 

17 Note (Xj_ j/g) reference (17J — 6 J _ 1 , (kj.}.^) reference [17] 

= X, 
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given contact angle This program also calculates contact point 
constant T but limits the contact angle to possibly 2 deg or 
greater 

Calculation of Natural Frequencies, Slosh Masses, and Their Location 
The remaining steps are relatively routine and, therefore, will not 
be described It is, however, lemarked that trapezoidal rule and 
midpoint rules weie employed conveniently in evaluating the in- 
tegrals For some quantities, quadratic fittings were mcde, 
such as dF/dR, before enteung the quadiature foimulas 


subject to 


#9 


+ Ah/ = 0, 

- 2m - 1 

X„ 2 W 

(40) 

II 

O 

at a: = 0 

(41) 

i- r ' 

at x = I 

(42) 


We shall select a complete set in (0, 1) which satisfies y = 0 at 
x = 0 but y 9^ 0 at x = 1 To be specific, we expiexs 


APPENDIX B 

A Modified Gctlerkin Method 

In the Galeikin method, it is generally assumed (191 that each 
of the cooidmate functions of a complete set satisfies the same 
boundary conditions as the exact solution This condition seems 
to be a sufficient condition lathei than a necessary condition 
since the series expansion as a whole may satisfy the prescubed 
conditions, not necessarily term by term Since the solution, if 
continuous, can be expanded into a conveigent series containing 
the complete set, the failure of the error minimization process, 
if it occurs, must he in the insufficient differentiability of the 
series Therefore, the present method as applied to a second- 
order ordinary differential equation performs integration by part, 
then imposes the remaining boundary condition 18 as illustrated 
by the following example 

Let y be governed by the simple equation 

18 This may be a new technique to satisfy the boundary condition 
or conditions However, a general proof and extension remain to 
be done 


where 


y = £ c « " m (*-?> 

m = 1 


„ 1 2»i — I 

A„ = : — it 


(43) 


Then, the modified Galerkin proceduie minimizes the error of the 
differential equation as follows 

0 - J sin (A,z) + X*yJ dx 

= ^sm (X„z) - j X„ ^ cos (\„x)dx 

+ J Xh/sin (\„x)dx = (— l)” -1 ?! 
nl eo 

- I X) C0S C0 * 

' JO m= 1 





n 1 oo 

H- I X 2 V) C m sin (X n a;) &in (K n x)dx 

JO T7I = 1 

*"Ti - (y - y) C n , n = 1, 2, , - (44) 


2 r,(- i r i {45) 

n (x„»-xn ( > 

The exact solution of the problem is 

y = r— — — r sm (\x) = C„' sin (\„x) (46) 

X cos X "i 

It is easy to show that the Fourier coefficient C„' of the exact 
solution is equal to C„, thus the modified Galerkm method yields 
the exact solution m the limit 

The late of convergence of the foregoing example (I\ 0) is 

oulj 1/n’ near the end, x = 1 This could be expected as dy/dx 
is /eio term by term while with infinite terms converge to T, as 
a whole when x — ► 1 Equation (13) in the text is analogously 
deuved, however, the rate of convergence might be better as illus- 
tiated by results in Appendix III. It should be pointed out that 
at the high Bond numbers, each auxiliary characteristic function 
is a noi mal mode and, as Bond number decreases, more character- 
istic functions are required to represent each mode Only by 
actual computation is one to find out how many of them are suf- 
ficient foi engineering purposes 


APPENDIX C 

Constant Contact Angle Condition 

For small perturbed motions at the instantaneous contact 
point, Part 2, Fig 5. 

dz = ft + /3 (47) 

= B ° + (S)„ dSr + On (48) 

df dr 6i / dr\ 

With tan do = — = — , 6i = sgn I — I, one has 

dr ds Vl + / r l \<y/ 

( a 4 /' ) -T (49) 

\ ds /II \dr l "\/l + /,V H L(! + /r ! ) /: Jn 



TANK WALL 


dsp - dsw cos 6c * hn cot ec 
. bn 
ds W" sfiTeF 

Fig 5 Geometry for contact pom! condition 


= ill'll + [" dsw 

L(1 + z r ») "Jb',11 


£t = * g “ (I) 


Similarly, at the instantaneous contact point 2, 


We shall assume constant contact angle independent of time l, 


6 C = dim — do = du, — d. = const (52) 


then, with dSiy = dS t — - dS,,- cos 8 C> Fig 6 

sm d„ 

— /rr£l ,, cos e c + ZRRt - „ = 0 atI1 

[1 -f /,*) <■ sm 8 C ds [1 + 2 fi ! ) ft sm 6 e 


i)h 

— = 7 h 
os 


1 T Z " e - a /rrft 1 

7 “ sin d« La + 2/) V '- C0S 9 ‘ (1 + /, s ) V, Jii 

For a convex tank with convex mtei face 

7 = ^ ~8 C ( (1 + n, _ C ° S Bc (1 +/, 2 )' /! lnj ( ° 6) 


8 "' = 6 " n + iW 


( 50 ) 
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APPENDIX III 
NUMERICAL EXAMPLES 
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Flat Interface with High Bond Number in a Cylindrical Tank 

N-g = 1000, —2- = 2. 34, a = 0. 68 with F = 0, a 12 X 18 mesh yielded 
to j a 


g 

dix B). 


= 1. 85 compared with 1. 847 from exact theory (Ref. 1 of Appen- 


m. 


•= 0. 193 compared with 0. 194 from high-G theory (Ref. 19 of 


m F 
Appendix B). 


B. 


D. 


= -0. 729 in. compared with -0. 724 in. from high-G theory (Ref. 19 
of Appendix B ). 

Flat Interface with Low Bond Number in a Cylindrical Tank = 10, 

h Q to? a 

= 2. 34, a = 0. 68 with T = 0, 12X18 mesh yielded = 2.15 com- 

g 

pared with 2. 46 from exact theory. A finer mesh is required for 
better agreement. 

h D 

Curved Interface with Low Bond Number N-o = 100, = 2. 34, 

2 H 7 a. 

^ , toia mi 

a = 0. 68, 9 = 2° yielded — =-= 1. 860, — ± = 0. 442, z. = -0. 724 in. 

g 2 P** 1 

oqa m| 

compared with theoretical value of = 1. 777, = 0. 438, 

g pa.3 

= -0. 73 in. (Ref. 12 of Appendix B ). It is noted that experimental 
to j a 

value of — — for 9 C - 0° is around 1. 78 to 1. 80. Finer net may lead 
to better approximation. 

k 

Spherical Tank with Flat Interface N-g = 1000, — = 1, a = 1 with T = 0, 

2 a 
CO -^3. 

11 X 11 mesh yielded — = 1. 54 compared with high-G sloshing value 

g 

to-, a 

of = 1. 54. 

g 
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E. 


*• Ai — ^ 

Spherical Tank with "Folding" Interface N R = 10, = 1. 14f, a = 1, 

co-i a mi 

e „ = 4. 96° (r = -65. 788), 12 X 12 mesh yield 1. 469, — 4 = 0. 604, 

c g 

Zi = -0.742 in. compared with Lockheed result (Ref. 16 of Appendix B) 
2 

CO jcl 

of = 1. 507 (for 6 = 5°). The first slosh mass is much larger than 

g c 

that of Lockheed, probably due to the inclusion of force directly 
attributed to surface tension in the present theory which yielded con- 
firmed results for cylindrical tanks (Case 3). The value of z^ would 
be the location of the center of the sphere if the contribution due to 
direct surface tension [ i. e. , the integral over F m Eq. (30b)] is 
neglected. 

Spheroidal tank with Folding Interface Ellipticity e = 0. 5, major axis 


a = 1, minor axis b = 0. 866, Ng =5, [3 = 0. 6263 (3/8 full tank), 
h n 

~= 0.4999, 6=5. 06°, (F = -38.430, V T = 1. 35)4 12 X 9 mesh 

3. C J-f 

2 

w l a tt m l 

(non-uniform on free surface) yielded =1.114 , = 0.582, 

g 2 PV L 

oo^a 

z-, = -0. 5 841 compared with Lockheed result of — = 0. 966. It is 

A & 

noted that the value of zj puts the first mass under the tank bottom 
which should be examined by future experiments when available. 


TCalculated for 3/4-full tank. 

JFor 0 C = 5.06° the Tused in ASME paper APM-EEE was in error and the 
correct value is given here as are the results. 

«? a m-, 

ttWith 23 X 17 mesh, = 1. 079, -r~- = 0.5736, z, = -0.4644; and 

a pV L 1 

2 2 

03^3. CO ■p 3 

= 15. 03 while Lockheed's result for is 13.26. It is uncertain which 

g g 

results are more accurate. The value of r has not been given by Lockheed 

to facilitate the explanation of the differences. 
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APPENDIX IV. 

LISTING AND INPUT SAMPLE OF TRIPOT 
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PROGRAM TRIPOT (INPUT, OUTPUT »TAPE7,TAPE6) 

C NORMAL COORDINATES 

C 02-1846-02 TRIANGULAR MESH GENERATOR AND SOLUTION OF 

C NORMAL MODES 

COMMON C (40 » 24 * 8) »PHI (40,24) , WB ,EPS , IT A , N , M , WO ,RHO ,R ( 6 ) ,YT(7) 
DIMENSION X (40,24) , Y (40 ,24) , ITITLE ( 12) 

DIMENSION ELAMDA (6) »SS(6),T(6),CT(6),CS(6),W(6),XT(7),A(6) 
DIMENSION CV (24) ,F (24,24) »EVAL (24,24) ,EVEC (24,24) ,EVI (1,1) 
DIMENSION EIGVAL (24) ,NVAL(24) *PHIS(24»24) »PHIW(64,24) 

DIMENSION DFDR (24) ,DRS(24) ,DFRR(24> ,DS(24) ,DRW{64) ,DZW(64) 
DIMENSION ALPJS (24) *ENU (24,24) ,GJ (24) »TKJ.(24) ,RMX (24,24) 

DIMENSION AMU (24) ,TALSK ( 24 ) »COMS ( 24) ,RW (64) ,ZW (64) ,EKS (24) ,TK(24) 
DIMENSION ELK (24) ,EMK(24> ,ZK(24) , DELI J (24,24) , I ROW (25) ,ICOL(25) 
DIMENSION SIGK (24) 

DIMENSION VNB ( 24 ) , VNF (24) ,VNW (40 ) ,XK(24) 

EQUIVALENCE (DFDR ( 1 ) ,C ( 1 ) ) , (DRS ( 1 ) »C (25) ) , (DFRR ( 1 ) ,C (49) ) , 

1 ( DS < I ) » C ( 73 ) > ♦ (DRW ( 1 ) , C (97) ) , (DZW (1 ) *C ( 161) ) , (ALPJS(l) ,C(225) ) , 

2 (ENU(l) ,C (249) ) ♦ (GJ(1) ,C(825) ) , (TKJ ( 1 ) , C (849) ) , (RMX ( I ) ,C ( 873) ) , 

3 ( AMU ( 1 ) , C ( 1449 ) ) , ( TALSK ( 1 ) , C ( 1473 ) ) , (COMS ( 1 ) ,C ( 1497) ) , 

4 (RW ( 1 ) , C ( 1521 ) ) , (ZW(1) »C( 1585) ) , (EKS(l) ,C(1649) ) , (TK (1) ,C(1673> ) , 

5 (ELK Cl) ,C (1697) ) , ( DELI J ( 1 ), C ( 1721 ) ) , ( IROW ( 1 ) , C ( 2297 ) ) , 

I 6 ( I COL (1) »C(2322) ) , (SIGK(l) ,C(2347) ) , (EMK<1) »C(2371) ) , 

! 7 ( ZK ( 1 ) , C ( 2395 ))»(XK(1)»C (2419 ) ) 

KUN=6 

LUN=7 

NDIM=24 

RHQ=0 ,05 

BETAC=0,5 

W0=.01 

Pl=3. 14159267 
READ 94, ITITLE 

94 FORMAT (12A6) 

READ 95,M,N» IT A , I RS , WA , WB , EPS , RLGTH , BOND 

95 FORMAT (4I5»5F10.0) 

PRINT 96, ITITLE, M,N, ITA,WA,WB,EPS, IRS 

96 FORMAT (26H1 TRIANGULAR MESH GENERATOR/, IX, 12A6/, 

117H LOGICAL MESH IS ,I5,3H X ,I5/,32H MAXIMUM NUMBER OF ITERATIONS 
2 IS, I5/»40H RELAXATION FOR LINEAR APPROXIMATION IS ,F7.3/, 

353H INITIAL RELAXATION FACTOR FOR NONLINEAR SOLUTION IS »F7.3/» 
428H EPSILON FOP CONVERGENCE IS ,F10.6/»19H RESTART IS NUMBER »I5/) 
NM1=N— 1 
MM 1 =M—1 
NM2=N-2 
EM=1 , 0 
MM=1 

NPMM1=N+M-1 

IF (MM ,EQ. 0) NM2=NM 1 
DO 97 J=2 , NMl 
DO 97 1=2, MM1 
C(I,J,1)=0.5 
C ( I , J*2) =0 ,5 
C ( I , J,3) =0 .5 

97 CONTINUE 

IF (IRS .EQ. 0) GO TO 3003 

GO TO (3001,1061,3006,3008,3012,377,389) , IRS 
3003 CONTINUE 

READ 200, (X(1»J) »Y (1 »J) , J=1 ,N) 

READ 200, (X(I,N) ,Y (I,N) ,1=2, MM 1) 

READ 200, (X(M,J),Y(M,J),J=1,N) 

READ 200, (X(I,1),Y(I,1),I=2,MM1) 

200 FORMAT (8F10.0) 

DO 93 J=1,N 
X(1,J)=X(1,J)/RLGTH 



Y<1,J)=Y<1»J)/RLGTH 
X (M*J)=X(M*J)/RLGTH 
93 Y(M*J)=Y(M»J)/RLGTH 
DO 92 I=2*MMl 
X(I»N)=X(I»N)/RLGTH 
Y(I,N)=Y(I»N)/RLGTH 
X(I,1)=X(I»1)/RLGTH 
92 Y(I,1)=Y<I*1)/RLGTH 

C INTERPOLATE FROM BOUNDARIES TO GET FIRST GUESS 

DO 100 J=2,NM1 
DO 100 1=2 »MM1 

X(I,J)=X(I,J-l)*X(I-liJ)-X(I-l,J-l) 

Y ( I » J) =Y { I » J— 1 )*Y(I-1»J)-Y(I-1 » J-l ) 

100 CONTINUE 
GO TO 3002 
C RESTART 1 

3001 READ (LUN) < (X < I * J) »Y ( I » J) , 1=1 ,M) * J=1 »N) 

WRITE (KUN) ( <X(I,J) ,Y(I»J) *I=1,M) *J=1,N) 

3002 CONTINUE 
IDIR=-1 

DO 103M0ST=1»ITA 

SX=0.0 

KIT=M0ST 

SY=0 • 0 

SRX=0 . 0 

SRY=0.0 

IDIR=-IDIR 

DO 102 K=2 *MMl 

DO 102 L=2 *NM1 

IF<IDIR) 1001*1002*1002 

1001 I=MM1-K*2 
J-NMl-L+2 
GO TO 1003 

1002 I=K 
J=L 

1003 CONTINUE 

XTEMP= (X ( I — 1 * J— 1 } +X ( I + l » J + l > +X ( I » J-l ) +X ( I * J+l ) +X < I — 1 » J) *X ( I + 1 * J) 

1 ) /6 • 0 

YTEMP=(Y{I-1*J-1)+Y(I + 1»J + 1)+Y(I»J-1)*Y(I»J+1)+Y(I-1*J)+Y(I+1»J) 

1 )/ 6.0 

RX=(X(I»J)-XTEMP)*WA 
RY=(Y(I»J) - YTEMP ) *WA 
X(I,J)=X(I,J)-RX 
Y ( I * J ) =Y ( I * J ) -RY 
SX=SX+X(I*J)**2 
SY=SY+Y(I»J)**2 
SRX=SRX+RX*»2 
SRY=SRY+RY**2 

102 CONTINUE 

IF (MOST-20* (MOST/20) ) 1022*1021*1022 

1021 WRITE (KUN) ( ( X ( I * J) ♦ Y ( I » J) * 1 = 1 *M ) , J=1 < N) 

REWIND KUN 

IRS=1 

IF(MOST .EQ. 20) PRINT 410, IRS 
410 FORMAT (21HOYOU MAY USE RESTART *14) 

1022 CONTINUE 
SRX=SQRT (SRX/SX) 

SRY=SQRT (SRY/SY) 

IF (SRX .LE. EPS .AND. SRY .LE. EPS) 105,103 

103 CONTINUE 

PRINT 201, SRX, SRY 

201 FORMAT (51H1PR0CESS DID NOT CONVERGE FOR INITIAL APPROXIMATION/, 

1 1 OH EPS-X IS »E15.7,12H EPS-Y IS ,E15.7//,41H X AND Y VALUES ARE 
2 PRINTED BELOW BY ROWS/) 

DO 104 1=1, M 

104 PRINT 202, (X (I,J) ,Y(I»J) *J=1*N) 
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202 FORMAT (8E15.7) 

WRITE (KUN) < <X(I,J) »Y(I*J),I=1»M) ,J=1,N> 

REWIND KUN 
IRS=1 

PRINT 410 » IRS 
STOP 

105 PRINT 203, KIT 

203 FORMAT (58H1PR0CESS CONVERGED FOR INITIAL APPROXIMATION ON ITERATI 
ION ,I5/»33H X AND Y ARE PRINTED BELOW BY ROW/) 

DO 106 1=1, M 

PRINT 202, { X ( I « J ) * Y ( I » J) »J=1»N) 

106 CONTINUE 

WRITE (KUN) ( (X(I,J) ,Y(I,J) ,I=1,M) ,J=1,N) 

1 » ( ( ( C ( I » J»K ) * 1 = 1 »M ) ,J=1*N) *K=1»3) 

REWIND KUN 
IRS=2 

PRINT 410, IRS 
GO TO 1062 
C RESTART 2 

1061 READ <LUN> ( (X < I , J) *Y ( I , J) , 1=1 ,M) , J=1 ,N) 

1 » ( ( (C(I*J»K> ,I=1,M) , J=1,N) *K=1 *3) 

WRITE (KUN) ((X(I,J),Y(I,J),I=1,M),J=1,N) 

1 » ( ( (C(I,J,K> »I=1,M) , J=1,N) *K=1 ,3) 

1062 CONTINUE 
WX=WB 
WY=WB 
IDIR=-1 

00 135 MOST=l , ITA 

IDIR=-IDIR 

K IT=M0ST 

SRX=0 • 0 

SRY=0.0 

SX=0.0 

SY=0 » 0 

DO 110 1=2, MM1 
DO 110 J=2,NM1 

DXDX=( <X(1-1*J-1)+2.0#X(I-1»J) +X(I,J+1) )-(X(I»J-l ) *2 . 0*X ( 1+1 » J) 

1 +X(I*1*J*1) ) )/6.0 

DYDX= ( (Y(I-1, J-l) +2. 0*Y ( 1-1 , J) +Y(I»J+1))-(Y(I»J-1) *2.0*Y ( I + 1 ♦ J ) 

1 +Y (I+1,J+1) ) )/6.0 
GAMMA=DXDX**2+DYDX**2 

DXDY= ( (X ( 1-1 , J) +2.0*X CI»J+1 > +X(I + 1,J+1))-(X(I-1 ♦ J-l ) +2 .0*X ( I ♦ J-l ) 

1 +X(I+1,J) ) )/6.0 

DYDY=( <Y(I-1*J)+2.0*Y(I,J+1 )*Y(I*1,J+1) ) - ( Y ( 1-1 ♦ J-l ) +2 .0*Y ( I , J-l ) 
1. + Y ( I + 1 , J ) ) )/6.0 
ALPHA=DXDY**2+DYDY**2 
BETA=DXDY«DXDX+DYDX*DYDY 
CPI =ALPHA-BETA 
CP2=BETA 
CP3=GAMMA-8ETA 

CPI =8E T AC^CPl ♦ (1.0-BETAC)*C(I,J,1) 

CP 2 = BE T AC*CP2+ ( 1 • 0-BET AC)*C(I,J,2) 

CP3=BE1 AC*CP3+ <1 .0-BETAC)*C(I , J,3) 

C ( I , J,1)=CP1 
C ( I , J,2) =CP2 
C ( I , J , 3 ) =CP3 
110 CONTINUE 

DO 120 K=2,MM1 
DO 120 L=2,NM1 
IF(IDIR) 1011,1012,1012 

1011 I=MMl-K+2 
J=NM 1 -L+2 
GO TO 1013 

1012 I=K 
J=L 

1013 CONTINUE 
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Cl=C(I,J,l) 

C2=CU,J,2> 

C3=C(I,J»3> 

C4=ci 

C5=C2 

C6=C3 

SUMC=C1 +C2+C3+C4+C5+C6 
IF (ABS(SUMC) - 1.0E-10) 120,120,111 
111 XTEMP=(C1*X(I-1»J) +C2*X ( 1-1 » J-l ) +C3#X ( I » J-l ) + C4#X { 1 + 1 » J) 

1 + C5*X(I + l,J+l)+C6*X(I,J-*-l) )/SUMC 
YTEMP=(C1*Y (I-l, J) +C2#Y (r-l*J-l> +C3*Y ( I * J-l ) + C4*Y ( I + 1 ♦ J ) 

1 +C5*Y(I+1,J+1)+C6*Y(I,J+1) )/SUMC 
RX= (X { I » J) -XTEMP) *WX 
RY=(Y( I»J)-YTEMP)*WY 
xti*j)=xa,j)-RX 
Y(I, J)=Y(I, J)-RY 
SRX=SRX+RX«*2 
SRY=SRY+RY**2 
SX=SX+X(I,J)**2 
SY=SY*Y(I,J)**2 
120 CONTINUE 

IF (MOST-20* (MOST/20) ) 3005,3004,3005 

3004 WRITE (KUN) ( ( X ( I , J ) , Y ( I , J) , 1=1 *M ) , J=1 ,N) 

1 , ( ( (C(I,J,K) ,I=1,M) »J=1,N) »K=1*3) 

REWIND KUN 
IRS=2 

IF(M0ST .EQ. 20) PRINT 410, IRS 

3005 CONTINUE 

IF (KIT-1) 126,126,127 

126 SPRX=SRX 
SPRY=SRY 

127 £PSCX=SQRT (SRX/SX) 

EPSCY=SQRT (SRY/SY) 

IF ( EPSCX .LE. EPS » AND. EPSCY .LE. EPS) GO TO 140 
ETAX=SQRT(SRX/SPRX) 

ETAY=SQRT(SRY/SPRY) 

ELX= (WX+ETAX-1 .0)/ (WX*SGRT (ETAX) ) 

ELY= (WY+ET AY-1.0 )/(WY*SQRT (ETAY) ) 

IF (ABS(ELX)-l.O) 129,129,128 

128 WAX=WX 

GO TO 130 

129 WAX = 2. 0/(1 . 0 + SQRT < 1 • 0-ELX#*2) ) -WO 

130 IF (ABS(ELY)-l.O) 132,132,131 

131 w A y = WY 

GO TO 133 

132 wAY=2.0/< 1 .0+SQRT ( 1 .0-ELY**2) ) -WO 

133 WX=RHO--WAX+(1.0-RHO)*WX 
WY =RHO*W AY + ( 1 .0~RhiO)*WY 
SPRX=SRX 

SPRY = SR Y 

135 CONTINUE 

PRINT 204, KIT, EPSCX, EPSCY 

?0<- FORMAT (43H1N0NL INEAR SOLUTION DID NOT CONVERGE AFTER ,I5,11H ITER 

1 * i IONS/, 1 OH EPS-X IS *E15.7 , 12H EPS-Y IS ,E15.7/,34H X AND Y ARE 

2 PRINTED BELOW BY ROWS/) 

DO 136 1=1, M 

PRINT 202, (X(I,J) »Y (I,J) *J=1,N) 

136 CONTINUE 

WRITE (KUN) ( (X(I,J) ,Y(I,J) ,1=1, M) ,J=1,N) 

1 , ( ( <C(I»J»K)»T=1,M) ,J=1,N) »K=1»3) 

IRS = 2 

PRINT 410, IRS 
STOP 

140 PRINT 205, KIT 

205 FORMAT (42H1 NONLINEAR SOLUTION CONVERGED ON ITERATION, 15/, 

1 34H X AND Y ARE PRINTED BELOW BY ROWS/) 
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3006 

3007 

300 

301 


DO 141 1=1, M 

PRINT 202, (X ( I * J ) »Y { I , J) *J=1 »N) 

141 CONTINUE 

WRITE (KUN) ( (X(I, J) ,Y (I, J) .1=1 ,M) »J=1»N) 
1 i (( (C(I,J,K),I=1,M) » J=1 ,N) »K=1»3) 

IRS=3 

PRINT 410, IRS 
GO TO 3007 
RESTART 3 
READ (LUN) 

♦ ( ( <C(I,J 
WRITE (KUN) 

, ( < <C(I,J 
CONTINUE 
DO 301 J=1 
DO 301 1=1 
DO 300 K= 1 
C(I, J,K)=1 
CONTINUE 
C(I,J,7)=0 
CONTINUE 
DO 302 J=1 
C(1,J,1)=0 
C ( 1 , J,2) =0 
C ( 1 , J, 6) =0 
C ( M , J, 3) =0 
C (M, J,4)=0 

302 C(M,J,5)=0 
DO 303 1=1 
C(I»1*1)=0 
C ( I , 1 ,2)=0 
C ( 1 , 1 ,3) =0 
C ( I , N , 4) =0 
C ( I , N,b ) =0 

303 C ( I , N , 6) =0 
DO 350 1=1 
DO 350 J=1 
DO 305 K= 1 
A (K) =0.0 
R (K) =0.0 
SS { K ) = 1 • 0 
T { K ) =1 . 0 
CT ( K ) =0 • 0 
CS (K ) =0 . 0 

ELAMDA(K)=C(I,J,K) 

XT ( K ) =0 • 0 

305 YT (K ) =0.0 
X IEMP=A < I , J) 

YTEMP=Y ( I , J) 

IF ([-!) 307,307,306 

306 Xl { t )=X ( 1-1 ♦ J) 

Y f « 1 ) =Y ( 1-1 , J) 

307 IF U-iu 308,309,309 
30d X l (4 j=a( 1 + 1 , J) 

YI (4) =Y ( 1+1 , J) 

309 IF (J-l) 311,311,310 

310 XT (3)=X(I, J-l) 


( (X(I,J) ,Y(I,J) ,I = 1,M) ,J=1,N) 
K) , 1=1 »M) ,U=1,N) ,K=1 ,3} 

( <X(I,J) ,Y(I,J) ,I=1,M> ,J=1,N) 
K) * 1 = 1 ♦ M) , J=1 ,N) ,K=1 »3) 


YT (3)=Y(I, J-l) 

311 IF (j-N) 312,313,313 

312 XT (6 ) =A ( I » J+ 1 ) 

YT ( 6 ) =Y (I, J+l) 

313 IF (1-1) 316,316,314 

314 IF (J-l) 316,316,315 

315 XT ( 2 ) =X ( 1-1 ♦ J-l ) 

YT (2 ) =Y ( 1-1 , J-l ) 

316 IF ( I-M) 317,319,319 
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317 IF (J-N) 318*319, 319 

318 XT (5 ) =X ( I + 1 , J+ 1 ) 

YT (5>=Y<I+l*J+l) 

319 CONTINUE 

XT ( 7 ) =XT ( 1 ) 

YT ( 7 ) =YT ( 1 ) 

DO 323 K=1 , 6 

320 SS <K > =SQRT ( ( XT < K ) -XTEMP ) **2+ ( YT <K ) -YTEMP ) **2) 

IF ( SS ( K ) .EQ. 0) SS<K)=1.0 

IF (K— 6) 321,322,322 

321 T(K)=SQRT( (XT(K^l)-XT(K) ) **2+ ( YT (K+ l ) -YT (K ) )**2) 

IF (T (K) .EQ. 0) T<K)=1.0 
R<K)=(XTEMP+XT(K)+XT(K+l) )/3.0 

GO TO 323 

322 T (K ) =SQRT ( (XT ( 6) -XT ( 1 ) ) ^*2+ ( YT (6) -YT ( 1 ) ) **2) 

IF (T (K) .EQ. 0) T(K)=1.0 
R<K)=(XTEMP+XT(6)+XT(l) )/3.0 

323 CONTINUE 

DO 340 K = l*6 
IF(K-I) 324,324,326 

324 IF (ELAMDA(6) > 325,331,325 

325 COSS=(T (6)**2+SS(6)**2-SS ( 1 ) **2) / < 2.0*T < 6 ) »SS (6) ) 

GO TO 328 

326 IF(ELAMDA<K-1> ) 327,331,327 

327 COSS=(T <K-1 ) *-*2 + SS <K-1 ) **2-SS (K) **2) / (2. 0*T (K-l ) *SS<K-1 ) ) 

328 CS<K)=1.0-C0SS**2 

IF (CS (K ) ) 329,329,330 

329 CS(K> =0.0001 

330 CS(K)=C0SS/SQRT(CS(K) ) 

331 IF (ELAMDA(K) > 332,340,332 

332 XA=0.5»(XT(K) ♦XTEMP) 

XB= (XT (K) +XT (K+l)+XTEMP)/3.0 
XC=0.5*(XT(K*1)+XTEMP) 

YA=0.5*<YT(K)+YTEMP> 

YB= ( YT (K)+YT (K+l)+YTEMP)/3.0 
YC=0.5*(YT(K+l)*YTEMP) 

AB=SQRT ( (XA-XB ) *«2+ (YA-YB) **2) 

8C=SQRT ( (YB-YC)**2+ (XB-XC) **2) 

AD=SQRT ( (XA-XTEMP) ( YA-YTEMP) *«2) 

BD=SQRT( (XB-XTEMP)**2+ < YB- YTEMP ) **2 ) 

CD=SQRT ( <XC-XTEMP)**2+ (YC- YTEMP >**2) 

ARA=0.5* (8C*CD+BD) 

ARA=SGRT(ARA*(ARA-BC)*<ARA-CD)*(ARA-BD) ) 

ARB=0 .5* ( AB+BD+AD ) 

ARB^SQRT { ARB* < ARB-AB ) * ( ARB-BD ) * ( ARB-AD ) ) 

A(K)=ARA,ARB 
IF (K-6) 334,333,333 

333 COST= (T (6) **2+SS ( 1 ) **2-SS ( 6 ) **2) / ( 2 . 0«T (6 ) *SS ( 1 ) ) 

GO TO 335 

334 COST=(T<K>*-*2 + SS<K+l)**2-SS<K)**2>/<2.0*T (K)*SS<K+1) ) 

335 CT(K)=1.0-COST»«2 

IF ( CT (K ) ) 336,336,337 

336 CT(K)=0.0001 

337 CT ( K ) = COST/SQRT (CT (K) ) 

340 CONTINUE 
K=1 

IF (1-1) 341,341,342 

341 IF ( J .EQ. 1) SS { 3 ) =0 . 0 
IF ( J .EQ. N) SS (6 ) =0 . 0 
CV(J)=(.5*(SS(3)+SS(6) ) )*XTEMP 

342 DO 343 K=2,6 

W(K)=ELAMDA(K)*R(K)*CT<K)+£LAMDA(K-1)*R<K-1>*CS(K> 

343 W(K)=0.5*W(K) 

W ( 1 ) =ELAMDA < 1 ) *R ( l ) *CT ( 1 ) +ELAMDA (6 ) *R (6) *CS ( l ) 
WU)=0.5*W(l) 

IF (XTEMP) 345,344,345 
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344 EMR=0 • 0 
GO TO 346 

345 EMR=EM*EM/XTEMP 

346 SUMC=0.0 

DO 347 K=1 , 6 
C(I»J,K)=W(K) 

347 SUMC=SUMC+W(K) +A(k)*EMR 
C ( I * J»8) =SUMC 

IF (J-I) 348.348,350 

348 IF (MM) 349*350.349 

349 C<I,J,1>=0.0 
CII,J,2)=0.0 
C(I,J,3)=0.0 
C<I»J»4)=0.O 
C(I,J»5)=0.0 
C(I» J»6)=0.0 
C ( I » J*8) =1 .0 

350 CONTINUE 

C SET COEFFICIENTS FOR CONTACT POINT 

C(1»N*1)=0.0 
C(1,N,2>=0.0 
C(1,N»5)=0.0 
C(1,N»6)=0.0 

READ 200*GFL6»GAMMA»THETAC 
THETAC=.0174532925»THETAC 
R ( 1 ) =X ( 1 »N-2 ) 

R ( 2 ) =X ( 1 ♦ N-l ) 

R(3)=X(1»N) 

YTU)=Y(l»N-2> 

YT ( 2 ) =Y ( 1 »N-1 ) 

YT(3)=Y(l.N) 

CALL SFIT(FR.FRR.XTEMP.YTEMP) 

C IF GFLG IS 0 COMPUTE GAMMA 

IF (GFLG) 3508,3503.3508 
3503 R(3)=X(1.N) 

R(2)=X(2,N) 

R(1)=X (3,N) 

YT ( 3 ) -Y ( 1 , N ) 

YT (2 ) =Y ( 2, N ) 

YT ( 1 )=Y (3.N) 

CALL SPIT (ZR,ZRR,XTEMP,YTEMP) 

ZTEMP=COS (THET AC ) 

IF(ABS(ZTEMP).GE.1.0) GO TO 3506 
GAMMA=(ABS(ZRR)-ABS(FRR)*ZTEMP)/SIN<THETAC> 

GO TO 3508 
3506 GAMMA=0.0 

3508 CONTINUE 

3509 CONTINUE 

PRINT 4000 ♦ GAMMA 
4000 FORMAT (7H0GAMMA »El5.8) 

WRITE (KUN) GAMMA, (((C(I*J,K)»I = 1,M),J=1.N) ,K=1 ,8) , (CV ( J) , J = 1 ,N> 
11 = 1 
IRS=4 

PRINT 410, IRS 
NM1=N 
GO TO 3011 
C RESTART * 

3008 11=1 
NM1 = N 

READ (LUN) ( (X(I,J) ,Y(I.J) , 1 = 1 ,M) ,J=1,N) 

I , ( ( (C(I,J.K) * 1 = 1 *M) ,J=1*N) ,K = 1 *3) 

WRITE (KUN) ( (X(I,J) ,Y(I,J) ,1=1, ft) » J=I ,N) 

1 , ( ( (C(I,J,K) * 1 = 1 *M) ,J=1*N) »K=1 *3) 

READ (LUN) GAMMA, (((C(I»J»K)»I=1»M),J=1,N)»K=1,8).(CV(J) »J=1»N) 
WRITE (KUN) GAMMA, {((C(I»J»K>»I=1»M),J=1»N)»K=1,8),(CV(J),J=1»N) 

3009 READ (LUN) (F ( J, I I ) , J=1 ,NM1 ) 
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IF ( EOF »LUN ) 3011,3010 

3010 11=11+1 

WRITE (KUN) (F(J»II)*J=1»NM1) 

GO TO 3009 

3011 CONTINUE 

DO 355 J=1,N 
DO 355 1=1, M 

355 PHI ( I , J) =0 ,0 

IF (MM ,NE. 0) CV(1)=0.0 
DO 360 ICOUNT=II»NMl 
C ( 1 ♦ ICOUNT ,7) =CV ( ICOUNT) 

CALL SOREL(l) 

DO 356 J=1 , NM1 
F(J,IC0UNT)=PHI (1,J) 

356 CONTINUE 

WRITE (KUN) (F ( J, ICOUNT) , J=l,NMl ) 

CC1* ICOUNT, 7)=0.0 

360 CONTINUE 
IRS=5 
NM1=N 

PRINT 410*1 RS 
GO TO 3013 
C RESTART 5 

3012 READ (LUN) ( (X ( I , J) ,Y ( I * J) * 1=1 ,M> * J=1 *N) 

1 , ( ( (C(I*J*K) * 1 = 1 *M) ,J=1,N) »K=1*3) 

NM1 = N 

WRITE (KUN) ( (X(I,J) ,Y (I, J) ,I=1,M) »J=1»N) 

1 * ( ( (C(I,J,K) *1-1, M) »U=1*N) *K=1*3) 

READ (LUN) GAMMA, ( ( (C (I, J»K> ,1=1 »M> » J=1»N) »K=1 ,8) , (CV <J> * J=1 *N) 
WRITE (KUN) GAMMA, { ( (C(I,J*K) ,I = 1»M) ,J=1*N) ,K=1,8) , (CV(J) ,J=1,N) 

DO 3014 1=1, NM1 

3014 READ (LUN) ( F ( J* I > , J=1 *NM1 ) 

WRITE (KUN) (F ( J, I > » J=1 »NM1 ) 

3013 CONTINUE 

DO 361 1=1, NM1 
DO 361 J=1 *NM1 
EVAL(I*J)=F(I» J> 

361 CONTINUE 

CALL EIGEN(NM1»EVAL*£VI»£VEC»50,3,0,1»NDIM) 

DO 372 1=1 *NM1 
NVAL(I)=I 

IF (EVAL ( I * I ) ) 370,370,371 

370 EIGVAL( I)=10.0E+60 
GO TO 372 

371 EIGVAL(I)=1.0/EVAL(I,I) 

372 CONTINUE 

DO 375 1=1, NM1 
K=I 

DO 374 J=K»NM1 

IF (EIGVAL(I)-EIGVAL(J) ) 374,374,373 

373 XTEMP=EIGVAL(J) 

II=NVAL<J) 

EIGVAL(J)=EIGVAL(I) 

NVAL ( J) -NVAL ( I ) 

E IGVAL ( I ) =XTEMP 
NVAL(I)=II 

374 CONTINUE ' 

375 CONTINUE 

PRINT 400 , ITITLE 

400 FORMAT (1H1,12A6/) 

PRINT 401 

401 FORMAT (77H EIGENVALUES AND EIGENVECTORS ASSOCIATED WITH INFLUENCE 
1 COEFFICIENT SOLUTIONS) 

NM2=N 

DO 376 1=1, NM2 
PRINT 402,1, EIGVAL(I) 
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402 FORMAT '(6H0MODE »I3,14H EIGENVALUE »E15.8) 

K=NVALU) 

PRINT 202 » (EVEC ( J,K) » J=1 »NMl ) 

376 CONTINUE 

WRITE (KUN) (ElGVAL(I) *I-1*NM1) » (NVAL(I) » 1=1 »NM1 } , 

1 ( (EVECC I,J) ,I=l,NMl) *J=1»NM1) * (NVAL(J) ,J=l,NMl) 

11 = 1 
IRS=6 

PRINT 410, IRS 
NM1=N 
NM2=N-1 
GO TO 382 
C RESTART 6 

377 READ (LUN) {(X(I»J)»Y(I*J)»I=1,M)*J=1»N) 

1 » ( < (C(I*J»K) » 1=1, M) *J=1»N) »K=1,3) 

NM1=N 

NM2=N-1 

WRITE (KUN) ( (X(I,J),Y(I,J) , I=1,M) *J=1,N> 

1 » ( ( (C(I»J»K) * 1=1 »M ) »J=1»N> »K=1»3) 

READ (LUN) GAMMA, (( <C(I,J»K) *I=1»M) »J=1*N) ,K=1»8) , <CV(J) ,J=1*N) 
WRITE (KUN) GAMMA, (((C(I,J,K),I=1,M) ,J=1,N) ,K=1,8) , (CV (J) , J=I,N) 
DO 378 1=1, NM1 
READ (LUN) (F < J, I ) ,J=1,NMI) 

378 WRITE (KUN) <F.( J, I ) , J=1 ,NM1 ) 

READ (LUN) (EIGVAL(I)*I=1,NM1),(NVAL(I)»I=1,NM1>, 

1 ( (EVEC(I,J) ,I=l*NMl) ,J=1,NM1) , (NVAL(J) ,J=l,NMl) 

WRITE (KUN) (ElGVAL(I) ,1=1, NM1) , (NVAL(I) *1=1, NM1) , 

1 ( (EVEC(I,J) ,1=1, NM1) ,J=1,NM1) » (NVAL(J) ,J=1,NMI) 

379 READ (LUN)(PHIS(I*II)»I=1*N)*(PHIW(I,II),I=1,NPMM1) 

IF (EOFtLUN) 382,381 

381 WRITE (KUN) (PHIS (I, II) »I = 1»N) * (PHIW ( 1 , 1 1 ) * 1=1 »NPMM1 ) 

11=11+1 

GO TO 379 

382 CONTINUE 

DO 383 J=1,N 
DO 383 1=2, M 

383 PHI ( I , J) =0 .0 

DO 388 ICOUNT=II ,NM2 
L=NVAL (ICOUNT) 

DO 384 1=1, NM1 

384 PHI (1,I)=EVEC(I,L) 

CALL SOREL (2) 

DO 385 J=1,N 

385 PHIS ( J , ICOUNT ) =PHI ( 1 , J) 

DO 386 J=1,N 

386 PHIW (J, ICOUNT) =PHI (M,J) 

DO 387 K=1,MM1 
KPN=K+N 

L=M-K 

387 PHIW(KPN,IC0UNT)=PHI(L,N) 

WRITE (KUN) (PHISU, ICOUNT), 1 = 1, N) , (PHIW ( I , ICOUNT) * 1 = 1 ,NPMMl ) 

388 CONTINUE 
NM1=N 
NM2=N-1 
GO TO 391 

C RESTART 7 

389 READ (LUN) ( ( X ( I , J ) , Y ( I , J ) , 1=1 , M) , J=1 , N ) 

1 , ( ( (C ( I » J,K) * 1 = 1 »M) »J=1»N) »K=1 *3) 

NM1=N 
NM2=N— 1 

WRITE (KUN) ( <X(I,J) ,Y(I,J) ,I=1,M) ,J=1»N> 

1 , ( ( (C(I, J,K) ,1=1 ,M) , J=1 ,N) ,K=1,3) 

READ (LUN) GAMMA, ( ( (C (I, J»K) ,1=1, M) , J=1*N> ,K=1,8) , (CV (J) »J=1»N> 
WRITE (KUN) GAMMA, ( ( (C ( I , J*K) , 1=1 »M) ,J=1,N) ,K=1,8) , <CV ( J) »J=1,N) 
DO 390 I=L,NM1 
READ (LUN) (F ( J, I") »U=1 ,NM1 ) 
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390 WRITE(KUN) IF ( J* I ) , J=1 *NM1 ) 

READ (LUN) (EIGVAL(I) »I=1»NM1) » (NVAL(I) »I=1»NM1) , 

1 ( ( EVEC ( I » J ) * 1 = 1 »NM 1 ) »J=l»NMl) » (NVAL(J) »J=l*NMl> 

WRITE (KUN) (EIGVAL(I) »I=1*NM1> » (NVAL(I) »I=1*NM1> ♦ 

1 ( <EVEC(I,J) *1 = 1 *NMl) »J=1*NM1> * (NVAL<J) * J=l»NMl> 

DO 392 J=1,NM2 

READ (LUN) <PHIS(I*J) *I=1*N) * <PHIW{I,J) »I=1*NPMM1> 

392 WRITE (KUN) (PHIS ( I ♦ J) ♦ 1=1 *N) * (PHIW ( I , J) , 1=1 , NPMM1 ) 

391 CONTINUE 
IRS=7 

PRINT 410* IRS 

NM1=N-1 

NM2=N-2 

C READ OTHER INPUTS 

READ 4002,ZCG*VOL,DEN»RHOU 

4002 FORMAT (5E15.8) 

PRINT 400,ITITLE 

PRINT 4003»BOND*ZCG*VOL*DEN»RHOU 

4003 FORMAT (13HOBOND NUMBER ,E15.8* 7H ZCG *E15.8* 9H VOLUME »E15.8 

1 //*22H LOWER FLUID DENSITY *£15,8 *25H DENSITY OF ULLAGE FLUID 

2 »E15.8) 

ZCG=ZCG/RL6TH 

V0L=V0L/RLGTH**3 

MM1=M-1 

NM1=N-1 

DO 4110 J=2,NM1 
XTEMP=X (1*J+1)-X(1*J) 

IF (XTEMP.EQ.O.) GO TO 4108 
DZDR=(Y(1*J+1)-Y(1*J) )/XTEMP 
IF <DZDR) 41 04*41 02 ,41 04 
4102 SIGN=1.0 

GO TO 4106 

4104 SIGN=DZDR/ABSF (DZDR) 

4106 DZDNF=SIGN/SQRT ( 1 . 0+DZDR**2 > 

DRDNF=-DZDR*DZDNF 
GO TO 4109 

4108 DZDNF=0.0 
DRDNF=-1 *0 

4109 VNF ( J) =Y ( 1 » J) *DRDNF-X ( 1 ♦ J) *DZDNF 

4110 CONTINUE 
VNF ( 1 ) =0 . 0 

VNF (N) =Y ( 1 ,N> *DRDNF-X ( 1 ,N ) *DZDNF 

MP1=M+1 

MP2=M+2 

DO 4120 I=2*M 

XTEMP=X (MP2-I tN)-X (MP1-I*N) 

IF (XTEMP.EQ.O. ) GO TO 4118 
DZDR= (Y (MP2-I »N)-Y (MP1-I»N) )/XTEMP 
IF <DZDR)4114*4ll2»4114 
4112 S IGN=— 1 • 0 
GO TO 4116 

4114 SIGN=DZDR/ABS(DZDR) 

4116 DZDNW=-S IGN/SQRT ( 1 ,0+DZDR**2) 

DRDNW=-DZDR*DZDNW 
GO TO 4119 

4118 DZDNW=0.0 
DRDNW=1.0 

4119 VNW (MP2-I )=Y (MP2-I *N)*DRDNW-X (MP2-I *N) *DZDNW 

4120 CONTINUE 

VNW ( 1 ) =Y ( 1 , N) *DRDNW-X ( 1 , N ) *DZONW 
DO 4130 J=1,NM1 
XTEMP=X(M»J+1)-X(M»J) 

IF (XTEMP.EQ.O. ) GO TO 4128 
DZDR=(Y (M*J+1)-Y ( M » J) )/XTEMP 
IF (DZDR) 4124*4122,4124 
4122 SlGN=l .0 
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GO TO 4126 

4124 S IGN=DZDR/ABS ( DZDR ) 

4126 OZDNW=~SIGN*DZDR/SQRT ( 1 . 0+DZDR**2 ) 

DRDNW=-DZDR*DZDNW 
GO TO 4129 

4128 0ZDNW=0.0 
DRDNW=1.0 

4129 VNB ( J) =Y (M, J) #DRDNW-X CM»J)*DZDNW 

4130 CONTINUE 

VNB(N)=Y(M»N)*ORDNW-X (M,N)*DZDNW 
C ( 1 » 1 , 7 ) =0 . 0 

SS(3)=SQRT( <X(1,N)-X(1»N-1> ) *#2+ < Y < 1 , N> -Y (1 , N-l > )**2> 
SS ( 4 ) =SQRT ( (X(2»N)-X(1»N) )**2+ (Y(2»N)-Y<1»N)) **2 ) 
C<1,N,7)=.5*(VNF(N)+VNF(N-1) ) * .5*SS (3) *X ( 1 ,N) ♦ 

1 .5* ( VNW < 1 ) +VNW (2 ) ) ».5*SS (4 ) *X ( 1 , N) 

DO 4132 J=2,NM1 

SS ( 3 ) =SQRT { <X(1»J)-X(1»J-1) >*»2+<Y(l,J)-Y<l,J-l) )*»2) 
SS(6)=SQRT< <X(1»J,1)-X(1»J) )**2+(Y(l,J+l)-Y<l»J> )**2) 
C<l,J,7)=VNF(J)*.5*(SS<3)+SS(6n*X<l,J) 

4132 CONTINUE 

00 4134 I=2*MM1 

SS < 1)=SQRT< (X< I»N)-X(I-1»N) )**2+ <Y ( I,N)-Y (I-1*N) >**2> 
SS(4>=SQRT< (X(I*1*N>-X(I»N) ) **2* ( Y ( I+l *N) -Y { I *N) )**2) 
C(I,N»7)=VNWCI)*.5*(SSn)+SS(4) >*X(I*N) 

4134 CONTINUE 

SS { 1 ) =SQRT (<X(M*N)-*X(M-1*N)) **2 + ( Y ( M*N) -Y (M-l *N) >**2> 
SS(3)=SQRT< <X(M»N)-X(M»N-1) )**2+ <Y(M,N)-Y <M,N-1 ) )**2> 
C<M,N»7)=(.5*(VNW(M)+VNW(M-1> )*.5*SS<3) 

1 +.5*<VNB(N)+VNB(N-1>)*.5*SS(1) >*X(M»N) 

DO 4136 J=2,NM1 

SS ( 3 ) =SQRT ( <X(M» J)-X(M, J-l) ) *#2 + ( Y ( M» J ) -Y <M*J-1) )**2> 
SS ( 6 ) =SQRT ( (X(M»J+1)-X(M»J) )**2+ (Y(M*J+1)-Y(M*J) )**2) 
C(M,J,7)=VNB<J)».5*(SS(3)+SS(6) )*X(M»J) 

4136 CONTINUE 

SS (6 ) =SQRT ( (X(M»2)-X(M»1) ) **2 + (Y (M*2)-Y (M*l ) )**2) 
C(M»1»7 )=VnB<1)*.5*SS(6)*.25*<3.*X(M,1>+X(M»2) ) 

DO 4138 1=1, M 
DO 4138 J=1,N 
PHI (I, J)=0.0 
4138 CONTINUE 

CALL SOREL(l) 

XTEMP=0.0 
00 4140 J=1,NM1 

XTEMP=XTEMP*«5*CPHI<1* J)*VNF<J)*X(1,J) 

1 +PHI(1,J+1) *VNF ( J+ 1 ) tt X ( 1 , J+l ) )* 

2 SORT ( (Y ( 1 » J+l )~Y( 1 » J) )«-*2+(X(l,J+l)-X(l»J) )**2> 
4140 CONTINUE 

00 4142 1=1, MM1 

XTEMP=XTEMP+ .5* (PHI ( I »N) *VNW ( I ) *X ( I » N> 

1 +PHI (I+1,N)*VNW(I*1)*X(I+1*N) >* 

2 SORT ( (Y(I+1,N)-Y(I»N) ) **2 + (X ( I+l »N) -X ( I »N) )**2) 
4142 CONTINUE 

DO 4144 J=1 ,NM1 

XTEMP=XTEMP+.5*(PHI (M»J)*VNB<J)*X<M»J) 

1 +PHI (M»J+1)*VNB(J+1)*X(M,J*1> )* 

2 SORT ( <Y(M,J+1>-Y<M»J) > **2+ (X <M, J+l > -X <M« J) )**2) 
4144 CONTINUE 

SIF=XTEMP*PI/V0L 
DO 415 J=1,NM1 
DRS(J)=X(1,J+1)-X(1,J) 

415 CONTINUE 

DRS(N)=DRS(NM1> 

DFDR(1)=0.0 

DFRR(1)=2.0*(Y(1,2)-Y(1,1) )/DRS(l)**2 
DO 417 J=2*NM1 
R ( 1 ) =X ( 1 , J-l ) 
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R(2)=X(1»J) 

R ( 3) =X ( 1 » J+l ) 

YT C 1 ) =Y { 1 » J-l) 

YT (2) =Y (1 * J) 

YT ( 3 ) =Y ( 1 » J+l ) 

CALL SFIT(FR»YTEMP»ZTEMP*XTEMP) 

OFDR ( J)=XTEMP 
OFRR ( J ) =YT£MP 
XK ( J) =ZTEMP 

417 CONTINUE 
DFDR ( N ) =FR 
DFRR(N)=DFRR(N-1) 

XK ( 1 ) =XK (2 ) 

XK (N)=YTEMP 

00 418 J=1,NM1 

418 DS(J)=2.0*PI»(X<1, J+1)+X(1»J> )*.5 

1 *SQRT ( (Yll»J+l)-Y(l>J) )«*2+<X(l,J+l)-X(l»J) )*» 2 ) 

DS (N) =05 (NMl ) 

DO 419 J=1,NM1 
DRW(J)=X(M*J+1)-X(M,J) 

RW(J)=X(M*J) 

ZW{ J) =Y (M * J) 

419 DZW ( J ) =Y (M» J+ 1 ) -Y (M » J) 

RW(N) =X (N*N) 

ZW(N)=Y(M,N> 

DO 420 K=l,MMl 
ITEMP=M-K+1 

DRW(N+K-1)=X(M-K,N>-X<ITEMP»N) 

DZW(N+K-1)=Y(M-K,N)-Y(ITEMP»N) 

RW(N+K)=XIM-K»N) 

420 ZW ( N+K ) =Y (M— K»N) 

RW ( N+M) =X ( 1 *N) 

ZW(N+M)=Y(1»N> 

DRW (NPMM1 ) =DRW (N+M-2) 

OZW (NPMM1 ) = DZW (N+M-2) 

DRW (N+M) =DRW (NPMM1 ) 

DZW (N+M) =DZW (NPMM1 ) 

DO 422 J=1*NM2 
XTEMP=0 • 0 
DO 421 K=1,NM1 

421 XTEMP=XTEMP +0.5* (PHIS (K* J) **2+PHIS (K+ 1 , J) +*2 ) *DS (K) 

422 ALPJS(J)=XTEMP 
DO 423 J=1 »NM2 
DO 423 1=1, NM2 

ENU(I»J)=2.0*PI/ALPJS(I)*X(1*N)*PHIS(N* J >*PHIS(N,I)»EIGVAL(J) 

423 CONTINUE 

C READ NUMBER OF EIGENVALUES TO BE USED IN REMAINING CALCULATIONS 

READ 9S,NEV 
NM2=NEV 

DO 430 J=1,NM2 
DO 430 1=1, NM2 
EPS I J=0 • 0 
6ETAIJ=0.0 
DEL I J ( I » J ) =0 . 0 
GAMI J=0.0 
XTEMP=0.0 
YTEMP=0 ♦ 0 

IF(I.EQ.J) DELIJ( I , J)=l .0 
DO 428 K=l,NMl 

YPRP=(Y(1»K+1)-Y(1»K) )/DRS(K) 

EE=DRS (K) /ABS ( DRS ( K ) ) 

CURVA=1.0+YPRP**2 
CURVB=SQRT (CURVA) 

EPSIJ=EPSlJ+0.5*<PHIS<K*J)*PHIS<K,I)/< (X(1,K> +.00001) «*2> 

1 +PHIS(K+1, J)*PHIS (K+1,I>/(X(1»K+1)»*2) )*DS(K> 

DPSIIA=(PHIS(K+1, [ )-PHIS(K,I) )/DRS(K) 
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DPSIJA=(PHIS(K+1*J)-PHIS(K*J) >/DRS(K) 

ZTEMP=YPRP**2/(CURVA*(0.5*(X <1,K)+X(1,K+1) ) )*#2)+ 

1 0.5*(XK(K)*»2*XK(K+1)**2) 

XTEMP=XTEMP+DPSIIA*DPSIJA*0S (K)/CURVA 

YTEMP=YTEMP-.5*ZTEMP*(PHIS(K»I)*PHIS(K,J) +PHIS (K+ 1 , I ) »PHIS (K + 1 , J) ) 
1 *DS(K> 

BETAIJ=BETAIJ*0.5»(PHIS(K*I)*PHIS(K,J)+PHIS<K+1*I)«PHIS(K+If J) ) 

1 *DS (K)/CURVB«EE 

428 CONTINUE 

EPSIJ=EPSIJ*EIGVAL<J)/ALPJS<I> 

GAMIJ=(XTEMP+YTEMP)*EIGVAL(J)/ALPJS<I) 

BETAIJ=BETAIJ*EIGVAL(J)/ALPJS(I) 

430 RMX ( I , J)=-GAMMA#ENU (I , J) +GAMI J+EM*EM*EPSI J+BOND*BETAI J 

1 * ( PEN-RHOU ) /DEN 1 

CALL MATINV(DELIJ*IROW,ICOL»NM2,NDIM»1.0E-06) 

DO 432 1=1 »NM2 
DO 432 J=1*NH2 
XTEMP=0.O 
DO 431 K=1,NM2 

431 XTEMP=XTEMP+RMX(I,K)*DELIJ(K, J) 

432 EVEC ( I , J) =XTEMP 
DO 433 1=1 *NM2 
DO 433 J=1 »NM2 

433 RMX(I»J)=EVEC(I»J) 

KM2=NM2 * 

CALL EIGEN(NM2»RMX»EVI»EVEC,50,2,0,1,NDIM) 

C PUT EIGENVALUES AND VECTORS IN INCREASING ORDER 

DO 442 1=1 ,NM2 
NVAL ( I ) =1 

IF (RMX (1*1) ) 440,440*441 

440 KM2=KM2-1 
COMS(I)=10.0E+60 
GO TO 442 

441 COMS ( I ) =RMX (1*1) 

442 CONTINUE 

DO 445 1=1, NM2 
K=I 

DO 444 J=K * NM2 

IF (COMS ( I ) -COM3 ( J) ) 444,444,443 

443 XTEMP=COMS(J) 

II=NVAL ( J) 

COMS ( J) =COMS ( I ) 

NVAL(J)=NVAL( I) 

COMS ( I ) =XTEMP 
NVAL ( I ) — 1 1 

444 CONTINUE 

445 CONTINUE 

DO 447 J=1 ,NM2 
K=NVAL(J1 
DO 447 1=1, NM2 
F ( I , J) =EVEC ( I ,K) 

447 CONTINUE 
PRINT 4005 

4005 FORMAT (39H0 EIGENVALUES AND EIGENVECTORS AT LOW G) 

DO 448 1=1, NM2 
PRINT 402 , I , COMS ( I ) 

PRINT 202, (F(J,I) ,J=1,NM2) 

448 CONTINUE 
XTEMP=OFOR (N) 

RII=X(1»N> 

F I I=Y ( 1 , N) 

NPMM2=NPMM1-1 

L=N+M 

RWMAX=RW(1) 

ZWMAX=ZW(1) 

DO 4491 1=2, L 
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IF (RW(I)-RWMAX) 44 91,449,449 

449 RWMAX=RW < I ) 

ZWMAX=ZW(I ) 

4491 CONTINUE 
RFMAX=X(1»1) 

ZFMAX=Y(l*i) 

DO 4493 I'=2 »N 

IF (X ( 1 » I ) -RFMAX ) 4493*4492*4492 

4492 RFMAX=XU»I) 

ZFMAX=Y(1,I) 

4493 CONTINUE 

00 453 J=1,NM1 
XTEMP=0.0 

00 451 K= 1 * NPMM2 

IF (ZW(K)-ZWMAX) 450*450,4501 

450 SIGN=1.0 
60 TO 4502 

4501 SIGN=-1.0 

4502 IF (DRW ( K ) ) 4504*4503*4504 

4503 DZDR=1.0E+07 
GO TO 4505 

4504 QZDR=DZW ( K ) /DRW (K ) 

4505 YrEMP=SQRT(1.0+DZDR**2) 

IF (A8S(DZDR)~1.0E+05) 4506,4507,4507 

4506 DRDN=SIGN*OZOR/YTEMP 
GO TO 4508 

4507 DRDN=1.0 

4508 DZON=-SIGN/YTEMP 

XT ( 1 )=0.5« (RW (K) +RW(K+1) ) 

YT { 1 ) =0 ,5* ( ZW (K)+ZW(K+1) > 

XT (2)=0,5<MPHIW(K»J) +PHIW (K+l* J> ) 

YT (2 ) =SQRT (DRW (K) **2+DZW (K)**2) 

XTEMP=XTEMP+PI*XT ( I ) *XT (2) * ( YT ( 1 ) *DRDN-XT ( 1 ) *DZDN> *YT (2) 

451 CONTINUE 

C XTEMP NOW CONTAINS MUJ1 

00 452 K=1,NM1 

IF (OFOR ( K ) ) 451 1*4512*451 1 

4511 SIGN=DFDR{K)/A8S(0FDR(K) ) 

GO TO 4513 

4512 SIGN=1.0 

4513 lF(ABS(OFOR(K) )-1.0E+05> 4514,4515*4515 

4514 DRDNA=-SIGN*DFDR (Kl/SQRT < 1 • 0+DFDR (K>**2) 

GO TO 4516 

4515 DRDN A=-l . 0 

4516 OZONA=S IGN/SQRT < 1 . 0+DFDR (K) **2) 

IF (DFDR (K+l ) >4521,4522, 4521 

4521 SIGN1=DFDR (K + l ) /ABS ( DFDR ( K + 1 ) ) 

GO TO 4523 

4522 S I GN 1=1,0 

4523 IF (ABS (DFDR (K+l) )-l .0E+05) 4527*4528,4528 

4527 DRDNB=-SIGN1*DFDR (K+1)/SQRT<1 , 0+DFDR (K+l )**2) 

GO TO 4529 

4528 DRDNB=-1.0 

4529 DZDNB=SIGNL/SQRT (1 . 0+DFDR (K+l >**2) 
XT=SQRT(DRS(K)»*2+(Y(1*K+1)-Y(1*K) ) »*2 ) 

XTEMP=XTEMP+0.5*PI*(PHIS(K»J)*(Y(1,K)*DRDNA-X ( 1 » K ) *DZDNA ) *X ( 1 »K ) 

1 +PHIS(K+1*J)*(Y ( 1*K+1 )*DRDNB-X (1 »K+1 ) *DZDNB) *X ( 1 *K+1) )*XT 

452 CONTINUE 
AMU ( J) =XTEMP 

453 CONTINUE 

IF(NEV ,GT. KM2 ) NEV=KM2 

DO 461 K=1,NEV 

DO 461 1=1, NEV 

TA|_SK(K)=0.0 

DO 460 J= 1 , NM 1 

XTEMP=0.0 



YTEMP=0.0 

XT=0.0 

YT=0.0 

00 459 L=1 *NEV 

XTEMP=XTEMP-F(L*K)*(-EIGVAL<L)*PHIS(J*L)*SQRT<1.0+DFDR<J)»*2) ) 
YTEMP=YTEMP+F(L»I)»PHIS(J*L) 

XT=XT— F (L *K ) * (— EIGV AL (L) tt PHIS(J + l»L) *SQRT (1,0 +0FDR ( J+ 1 ) #*2 ) ) 

459 YT=YT+F(L»I)*PHIS(J*1*L> 

460 TALSK CK) = ( ,5*XTEMP*YTEMP»X U » J) +0 ,5*XT*YT*X ( 1 * J + l > ) *ABS (DRS ( J) ) 
1 +TALSK(K) 

461 RMX(I»K)a2.0*PI*TALSK(K) 

CALL MATINV (RMX » I ROW * ICOL*NEV »NDIM* 1 .OE-06) 

DO 464 K=1,NEV 

XTEMP=0.0 

DO 463 J=1 *NM1 

XTEMP=XTEMP + 0 .5* (-EIGVAL (K) *PHIS ( J »K) *X ( 1 » J) 

1 -EIGVAL(K)*PHIS(J+1»K>*X(1» J+l) )*DS(J) 

463 CONTINUE 

464 GJ(K)=-XTEMP 
DO 466 K=1,NEV 
TALSK(K)=0.0 
DO 465 J=1*NEV 

465 TALSK(K)=TALSK<K)+F<JfK)«GJ<J) 

466 CONTINUE 

DO 4670 1=1 *NEV 

XTEMP=0.0 

DO 467 J=1*NEV 

467 XTEMP=XTEMP+RMX ( I , J) »TALSK ( J) 

4670 EKS(I)=XTEMP 

EMS=1.0 

ZS=0.0 

DO 468 K=1*NEV 

468 SIGK (K) =GJ (K) 

DO 471 K=1»NEV 

XTEMP=0.0 

XT=0,0 

YT=0 . 0 

DO 470 J=1»N£V 
XTEMP=XTEMP+F ( J*K) *SIGK ( J) 

470 XT=XT+F(J»K)*AMU(J) 

TK (K)=0.5*XTEMP 
ELK (K > =EKS (K) *XT 
EMK(K)=EKS(K)*TK(K)/VOL 
EMS=EHS-EMK(K) 

ZK (K)=RLGTH/VOL^ELK (K)/EMK<K) 

ZS=ZS+EMK(K)*ZK(K) 

471 CONTINUE 
ZS=1»0/EMS* (ZCG-ZS) 

XTEMP=0 < 0 

DO 473 K= 1 » NE V 

473 XTEMP=XTEMP+EMK(K)*ZK(K)^*2 

SI0=SIF-(EMS*ZS**2+XTEMP)/RLGTH**2 

EMF=DEN*V0L*RLGTH**3 

PRINT 475»SIF*SlO*EMF 

475 FORMAT ( 25H0 SIF-MOMENT OF INERTIA =*E15.8/ 

1 25H0 S I O-MOMENT OF INERTIA =*E15.8/ 

2 22H0 EMF-MASS OF LIQUID =*E15.8) 

PRINT 500 

500 FORMAT ( 7H0 MK/MF/) 

PRINT 202* (EMK(K) ,K=1*NEV) 

PRINT 501 

501 FORMAT (4H0 ZK/) 

PRINT 202* (ZK(K) ,K=1*NEV) 

PRINT 502 * EMS » ZS 

502 FORMAT(5HOMO= *E15.8*7H Z0= *E15.8) 

ENDFILE KUN 
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STOP 

END 
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SUBROUTINE SFIT ( YR,XK3,XK,YR2) 

C FIT QUADRATIC THROUGH THREE POINTS 

COMMON C <40, 24,8 ), PHI (40*24) ,WB,EPS, ITA ,N»M, WO ,RHO,R ( 6 ) ,YT<7) 
IF(ABS(<YT<3)-YT(2))/<R<3)~R(2))).GT.1.0) GO TO 30 
IF(ABS( (YT(2)-YT(1) )/(R(2)-R(l) ) ) .GT.1.0) GO TO 30 
XTEMP=(R(1)**2-R(3>«*2)*(R(2)-R(3) )- (R ( 2) **2-R (3) **2) # <R (1 ) -R (3) ) 
YTEMP=(YT(1)-YT (3) ) * ( R (2 ) -R ( 3) ) - < YT (2) -YT <3> > * (R < 1 ) -R (3) ) 
AD=YTEMP/XTEMP 

AB=(YT(2)-YT(3)~AD*<R(2)**2-R(3)**2) )/(R(2)-R(3) ) 

YR=2.0*AD*R (3) +AB 
YR2=2.0*AD*R<2) +AB 
YRR=2.0*AD 

XK=YRR/ ( 1 .0+YR2**2) .5 
H2=R ( 2 ) -R ( 3 ) 

Hl=R ( 1 ) — R (3) 

F2=YT(2) 

FI=YT(1) 

F0=YT (3) 

YRR=2.*<H2*(F1-F0)— H1MF2— F0> ) / < H1*H2* ( H1-H2) ) 

XK3=YRR/(1.0+YR**2)**1.5 

RETURN 

30 YTEMP= ( YT ( 1 ) **2-YT (3) **2 ) * ( YT (2) -YT (3) ) - (YT (2) **2-YT (3) **2) * 

1 ( YT ( 1 ) — YT ( 3 ) ) 

XTEMP=(R(1)-R(3) )*<YT(2)-YT<3) )-(R(2)-R<3) >* (YT(l)-YT (3) ) 
AD=XTEMP/YTEMP 

AB=(R(2)-R(3)-AD»(YT(2)*»2-YT(3)**2) ) / ( YT (2) -YT ( 3) ) 
YR=1.0/(2.0*AD*YT<3) +AB) 

YR2=1 ,0/(2.0»AO*YT (2) +AB) 

YRR=-2.0*AD*YR2*»3 
XK=YRR/ ( 1.0+-YR2**2> **1 .5 
H2=YT(2)-YT (3) 

Hl=YT ( 1 ) -YT ( 3 ) 

F2=R(2) 

F1=R(1) 

F0=R(3) 

RYY=2.« <H2*CF1-F0)-H1*(F2-F0) ) / ( H1#H2* ( H1-H2) ) 

YRR3=-RYY*YR**3 
XK3=YRR3/ ( 1 •0+YR**2) **1.5 
RETURN 
END 



SUBROUTINE SOREL(LL) 

COMMON C<40,24,8) ,PHI(40,24> »WB,EP5» ITA»N»M»W0,RH0»R(6)»YT (7) 
SUCCESSIVE OVER-RELAXATION TO SOLVE FOR PHI 

il=ll 

WX=WB 

DO 135 MOST=l»lTA 

KlT=MOST 

SRX=0.0 

SX=0.0 

DO 120 I=IL*M 
DO 120 J=1 »N 
PHIT=C(I»J»7) 

DO 108 K=l,6 
CT=C ( I » J»K) 

IF (CT .EG. 0.0) GO TO 108 
IF (K-2 ) 100,101,102 
100 11=1-1 
JJ=J 

GO TO 107 
101 11=1-1 
JJ=J-1 
GO TO 107 

102 IF { K-4) 103,104,1 05 

103 1 1 = I 
JJ=J-1 

GO TO 107 

104 11=1+1 
JJ=J 

GO TO 107 

105 IF (K-5) 109,106,109 

106 11=1+1 
JJ=J+1 

GO TO 107 
109 II=I 
JJ=J+1 

107 PHIT=PHIT+CT*PHI (II, JJ) 

108 CONTINUE 
PHIT=PHIT/C(I»J»8) 

RX=WX*(PHI ( I , J ) — PHIT ) 

PHI ( I , J) =PHI (I,J)-RX 

SRX=SRX+RX**2 

SX=SX+PHI(I,J)»#2 

120 CONTINUE 

IF (SX .EQ. 0.0) GO TO 140 
IF (KIT-1) 126,126,127 

126 SPRX=SRX 

127 EPSCX=SQRT (SRX/SX) 

IF (EPSCX .LE. EPS) GO TO 140 
ETAX=SQRT (SRX/SPRX) 

ELX= (WX+ETAX-1 .0)/<WX*SGRT (ETAX) ) 

IF (ABS(ELX)-l.O) 129,129,128 

128 WAX=WX 

GO TO 130 

129 |*/AX=2.0/(1.0+SQRT(1.0-ELX**2) )-WO 

130 WX=RHO*WAX+<1.0-RHO)*WX 
SPRX=SRX 

135 CONTINUE 

IF(IL .EQ. 1) PRINT 200 , K IT , EPSCX 
IF(IL .EQ. 2) PRINT 204, K I T , EPSCX 

200 FORMAT <*1SOLUTION DID NOT CONVERGE FOR INFLUENCE COEFFICIENTS*/, 
1 * AFTER *,I5»* ITERATIONS EPSCX= *,E15.7) 

204 FORMAT (*1S0LUTI0N DID NOT CONVERGE FOR EIGENVECTOR ON SURFACE*/, 

1 * AFTER *,I5»* ITERATIONS EPSCX= *»E15.8) 
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DO 136 1=1, M 

PRINT 202, (PHI <I,J) ,J=1,N) 

202 FORMAT (8E15.7) 

136 CONTINUE 

STOP 

140 IF ( IL .EQ. 1) PRINT 203,KIT 
IF(IL .EQ. 2) PRINT 205iKIT 

203 FORMAT (* INFLUENCE COEF. SOL. CONVERGED AFTER ITERATION 1 *, 15) 
205 FORMAT (* EIGENVECTOR SOLUTION CONVERGED AFTER ITERATION 1 *, 15 ) 

RETURN 

END 
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SUBROUTINE EIGEN (N»B»TI *T» MAX IT »NDEC*N10PT » N2 OPT » NDIM ) 

SUBROUTINE FOR GENERATING THE EIGENVALUES AND EIGENVECTORS 
OF A REAL SYMMETRIC OR NON-SYMMETRIC MATRIX. 

THIS PROGRAM GENERATES THE EIGENVALUE MATRIX (REAL OR COMPLEX)* 
AND AS OPTIONS, THE EIGENVECTOR MATRIX AND ITS INVERSE. 

THE CALL FOR THIS SUBROUTINE IS AS FOLLOWS* 

CALL EIGEN (N*B*TI*T»MAXIT,NDEC»NI0PT»N20PT*NDIM) 

WHERE N IS THE ORDER OF THE MATRIX 

B IS THE MATRIX WHOSE EIGENVALUES ARE DESIRED 
TI IS THE INVERSE OF THE EIGENVECTOR MATRIX 
T IS THE EIGENVECTOR MATRIX 

MAXIT IS MAX NO. OF ITERATIONS FOR GENERATING EIGENVALUES 
NDEC IS THE NUMBER OF TIMES RESULT IS REFINED 
NIOPT IS 1 IF OPTION =1 IS DESIRED* OTHERWISE 0. 

OPTION =1 PROVIDES FOR GENERATING THE EIGENVECTOR MATRIX INVERSE 
N20PT IS 1 IF OPTION =2 IS DESIRED, OTHERWISE 0. 

OPTION =2 PROVIDES FOR GENERATING THE EIGENVECTOR MATRIX 
NDIM IS DIMENSIONED NO. OF ROWS OF MATRIX (B) 

THE ORGINAL MATRIX B IS LOST DURING THE COMPUTATIONS AND IS 
REPLACED BY THE EIGENVALUE MATRIX. 

DIMENSION B(1),TI(1),T(1) 

C INITIALIZE COUNTERS FOR NO. OF ITERATIONS AND YR» YS REDUCTIONS 

IT = 0 

NT IMES=0 
ANORM=0.0 
- DO 1100 1 = 1 *N 
DO 1100 J=1,N 
IJ=( J-1)*NDIM+I 

1100 ANORM=ANORM+B ( I J ) **2 
ANORM=SQRTF (ANORM) 

DO 1101 1=1, N 

DO 1101 J=1,N 
IJ=(J-1)*NDIM+I 

1101 8 ( I J) =B ( I J) /ANORM 

C FORM IDENTITY MATRIX IN TI LOCATION IF OPTION 1 IS DESIRED 

1000 IF(NIOPT) 1010*1010*1001 

1001 DO 1003 1=1, N 

II = < I — 1 ) *NDIM+ I 
DO 1002 J=1,N 
IJ = < J-l )*NDIM+I 

1002 TI (IJ) = 0. 

1003 TI ( II) = 1.0 

C FORM IDENTITY MATRIX IN T LOCATION IF OPTION 2 DESIRED 

1010 IF (N20PT) 1020, 1020*1011 

1011 DO 1013 1=1, N 

II = (I-1)#NDIM+I 
DO 1012 J=1,N 
IJ = ( J-l ) *ND IM+ I 

1012 T ( I J) = 0. 

1013 T ( 1 1 ) = 1.0 
1020 CONTINUE 

YR=10.0E-7 
YS=1 0 . OE-7 

96 NO=N-l 
SUM=1 0 . 0E20 

97 T AU=0 • 0 
EN=Q . 0 

DO 1 1=1, NO 
JO=I+l 
DO 1 J=JO» N 
IJ = (J-1)»NDIM+I 
JI = <I“1)*NDIM+J 
1 TAU = TAU+B ( IJ) **2+8 ( JI )*>2 
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DO 2 1 = 1 »N 

II = ( 1-1 ) *NDIM + I 

TE = B ( 1 1 ) 

2 EN=EN+TE*»2 

4 EN=EN+TAU 

5 DELN=SUM-EN 

6 IF (DELN) 8,8,7 

7 SUM=EN 
IT=IT+ 1 

IF (MAXIT-IT) 120,120,10 

8 CONTINUE 

IF { NDEC-NT IMES ) 120,120,9 

9 YR=YR/100.0 
YS=YS/100.0 
NTIMES=NTIMES+1 
IT=IT+1 

10 DO 98 K=1 » NO 

KK = (K-l ) *NDIM+K 
KO=K*l 

11 DO 98 M=KO,N 

MM = (M-1)*NDIM+M 

KM = (M-1)*NDIM+K 

MK = (K-1)*NDIM+M 

H=0 « 0 

G=0.0 

HJ=0.0 

12 DO 24 1=1, N 

13 IF(I-K) 14,24,14 

14 IF (I -M >15,24, 18 

15 IK = (K-1)*NDIM+I 
Kl = (I-1)*NDIM+K 
IM = <M-1)*NDIM+I 
MI = (I-1)*NDIM+M 
80 = B ( IK ) 

16 BR = B (KI ) 

17 8Q = BUM) 

18 6S = B (MI) 

19 H=H+BR*BS-BO*BQ 

20 TEP=B0*BO+BS*BS 

21 TEM=BR*BR+BQ*BQ 

22 G=G+TEP+TEM 

23 HJ=H J-l EP+TEM 

24 CONTINUE 

25 H=2.0*H 

D = B (KK) -B (MM) 

TEP = B (KM) 

TEM = B(MK) 

C=TEP+TEM 

E=TEP-TEM 

26 IF (ABSF (C)~YR)27,27,30 
27' CC=1.0 

28 SS=0.0 

29 GO TO 39 

30 8Y=D/C 

31 IF (BY)400»40 1,401 

400 SIG=— 1 .0 
GO TO 32 

401 SIG =1,0 

32 COT=BY+ (SIG*SQRTF(BY*BY+1 ,0) ) 

33 SS=SIG/SQRTF(COT*COT+1.0) 

34 CC=SS*C0T 

35 TEP=CC*CC-SS#SS 

36 TEM=2.0*SS*CC 

37 D=D*TEP*C*TEM 

38 H=H*TEP-HJ<UEM 

39 CONTINUE 
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40 ED=2.0*E*D 

41 EDH=ED-H 

42 DEN=G+2.0*(E*E+D*0> 

43 TEE=EDH/ (DEN+DEN) 

75 CONTINUE 

205 IF(ABSF(TEE)-YS>44,44,46 

44 CH=1.0 

45 SH=0.0 
GO TO 48 

46 CH=1.0/SQRTF(1.0-TEE*TEE) 

47 SH=TEE*CH 

48 C1=CH*CC-SH*SS 

49 C2=CH*CC+SH*SS 

50 S1=CH#SS+SH*CC 

51 S2=-CH*SS+SH*CC 

52 CONTINUE 

53 IF (Si ) 55*54*55 

54 IF(S2) 55*98*55 

55 DO 59 J=1»N 

KJ = (J-1)*NDIM+K 
MJ = (J-1)*NDIM+M 

56 BO = B (K J ) 

57 BR = B (MJ) 

58 B(KJ) = Cl*BO+Sl*BR 

59 B ( M J ) = S2'*80+C2*8R 

60 DO 66 J=1*N 

JK = (K-l ) *NDIM+ J 
JM = ( M— 1 ) *ND3M+ J 

61 BO = B ( JK) 

62 BR = B ( JM ) 

65 6 ( JK ) = B0*C2-BR*S2 

66 B ( JM ) = -80*S1+BR*C1 

1070 IF(NIOPT) 1075*1075*1071 

1071 DO 1072 J=1,N 

KJ = ( J-1)*NDIM+K 
MJ = ( J-l ) tt NDlM+M 
BQ = T I (KJ) 

BS = TI(MJ) 

T I (KJ) = C1*BQ+S1*BS 

1072 TI(MJ) = S2*BQ+C2*BS 

1075 IF(N20PT)98*98,1076 

1076 DO 1077 J=1»N 

JK = (K-l )«NDIM+J 
JM = (M-1)*NDIM+J 
BO = T (JK) 

BR = T ( JM ) 

T ( JK ) = BO*C2-BR*S2 

1077 T(JM) = -B0*S1+BR*C1 
98 CONTINUE 

GO TO 97 

120 DO 1102 1 = 1, N 
DO 1102 J=1»N 
IJ=(J-1)<*NDIM+I 

1102 B(IJ)=B(IJ)«ANORM 

IF ( N20PT ) 1107*1107,1103 

1103 DO 1106 J=1*N 
ANORM=0.0 

DO 1104 1=1, N 
IJ=(J-1)*NDIM+I 

1104 ANORM=ANORM + T ( IJ>*"*2 
ANORM=SQRTF (ANORM) 

DO 1105 1=1, N 
IJ=(J-1 )*NDIM+I 

1105 T < I J) =T' ( I J) / ANORM 

1106 CONTINUE 

1107 RETURN 
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SUBROUTINE MPRINT ( A*M»N*MD) 

DIMENSION A ( 1 ) » IT (6) »C (b) 

EQUIVALENCE (IT»C) 

2 FORMAT ( 1HO » AX, 6( 6X» 7HCOLUMN 114 ) 

3 FORMAT <1H 114, X, (6E 17.8) ) 

Nl=N 

N2=6 

N3=6 

N4=l 

A IF (N3-N1) 6,6,5 

5 N2=Nl-N3«-6 
N3=N1 

6 K-0 

DO 7 1= N4,N3 

K=K + 1 

7 IT ( K ) — I 

PRINT 2 » < IT ( I ) » 1 = 1 ,N2) 

DO 9 I-1,M 

K=0 

L=MD* (N4-1 ) +1 
DO 8 J=N4,N3 

K=K + 1 
C(K)=A(L) 

L=L+MD 

8 CONTINUE 

9 PRINT 3, I, (C(K) ,K=1»N2) 

IF (N3-N1) 10*11*11 

10 N3=N3+6 
N4=N4+6 
GO TO 4 

11 RETURN 
END 


/// ) 
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C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 



SUBROUTINE MAT INV (A* I ROW* ICOL, N»NDIM , SMLST ) 

MATINOOO 


DIMENSION A ( 1 ) * IROW ( 1 ) * ICOL ( 1 ) 

MATIN001 


709-16065 

MAT IN002 


709-16065 SUBROUTINE MAT INV - MATRIX INVERSION ROUTINE 

MATIN003 
MAT IN004 


A = ARRAY NAME OF MATRIX 

MATIN005 


IROW = DIMENSIONED AT N+l OR GREATER 

MATIN006 


ICOL = DIMENSIONED AT N OR GREATER 

MATIN007 


N = NUMBER OF EQUATIONS 

MAT IN008 


NOIM = VALUE OF |I IN DIMENSION A(I*J) * I AND J MAY DIFFER 

MATIN009 


SMLST = SMALLEST ‘LEADING ELEMENT ALLOWED BEFORE CALLING THE 

MATIN010 


SYSTEM SINGULAR * USUALLY =1.0 E-04 OR 1.0 E-05 

MAT I NO 1 1 
MAT INO 12 


NP1=N+1 

MATIN013 


DO 5 1=1, N 

MAT INO 14 


ICOL ( I ) -I 

MATIN015 

5 

IROW (I ) =1 

MATIN016 


DO 75 ITER=1,N 

MATIN017 


MAXR=ITER 

MATIN018 


MAXC=1 

MAT INO 1 9 


TEMP=ABSF(A(MAXR>) 

MATIN020 


LlMITC=NPl-ITER 

MATIN021 


DO 15 I=ITER,N 

MATIN022 


DO 15 J=1,LIMITC 

MATIN023 


IJ=(J-1)»NDIM+I 

MATIN024 

10 

IF ( TEMP- (ABSF ( A ( I j) ) ) ) 10,15,15 

MATIN025 

MAXR=I 

MATIN026 


MAXC=J 

MATIN027 


TEMP=ABSF(A(IJ) ) 

MAT IN028 

15 

CONTINUE 

MAT I N029 


IF (TEMP-SMLST) 20,20,25 

MATIN030 

20 

IROW (NPl ) =ITER 

MATIN031 


PRINT 200, ITER, SMLST 

MATIN032 

200 

FORMAT ( 7HQ0N THEI3,63HTH ITERATION ALL THE REMAINING TERMS WERE 

LMATIN033 

IESS THAN OR EQUAL TO E11.4,18H IN ABSOLUTE VALUE) 

MATIN034 


RETURN 

MAT IN035 

25 

IF (MAXR-ITER) 30,40,30 

MAT IN036 

30 

DO 35 J=1,N 

MAT IN037 


MAXRJ=(J-1)*NDIM+MAXR 

MATIN038 


ITJ=(J-1) *NDIM+ITER 

MATIN039 


TEMP=A (MAXRJ ) 

MAT IN040 


A(MAXRJ)=A(ITJ) 

MAT IN041 

35 

A ( IT J) =TEMP 

MATIN042 


ITEMP=IR0W (MAXR) 

MAT IN043 


IROW (MAXR)=IROW (ITER) 

MAT IN044 


IROW(ITER)=ITEMP 

MAT IN045 

40 

IF (MAXC-1) 45,55,45 

MATIN046 

45 

DO 50 1=1, N 

MAT IN047 


IMAXC= (MAXC-1 ) »NDIM+I 

MATIN048 


TEMP=A( I) 

MATIN049 


A ( I ) =A ( IM AXC) 

MAT IN050 

50 

A ( IMAXC) =TEMP 

MATIN051 


ITEMP=ICOL(MAXC) 

MAT ING52 


ICOL (MAXC ) =ICOI_ ( 1 ) 

MATIN053 


ICOL (1)=ITEMP 

MATIN054 

55 

TEMP=A(ITER) 

MATIN055 


ITEMP=ICOL(l) 

MATIN056 


DO 60 J=2 ,N 

MAT IN057 


IT JM 1= ( J-2 ) *NDIM + ITER 

MATIN058 


ITJ=(J-1)*NDIM+ITER 

MATIN059 


A(ITJM1)=A(ITJ)/TEMP 

MAT IN060 

60 

ICOL ( J— 1 > =ICOL ( J) 

MAT IN061 


ITN= (N-l )*NDIM+ITER 

MATIN062 
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A(ITN)=1. 0/TEMP 

MATIN063 


ICOL(N)=ITEMP 

MATIN064 


DO 75 1=1 *N 

MAT IN065 


IF ( I— ITER) 65,75,65 

MATIN066 

65 

TEMP=A(I) 

MATIN067 


DO 70 J=2,N 

MATIN068 


IJM1=(J-2)#NDIM*I 

MATIN069 


I J= ( J— 1 ) *NDIM+I 

MATIN070 


ITJM1=(J-2)*NDIM+ITER 

MATIN071 


A ( I JM1 ) =A ( I J) -A ( ITJM1 ) #TEMP 

MATINO 12 

70 

CONTINUE 

MATINQ73 


IN= t N — 1 ) *NDIM+ I 

MAT INO 74 


ITN=(N-1)*NDIM+ITER 

MATINO 75 


A ( IN > =- (TEMP^A ( ITN) ) 

MATIN076 

75 

CONTINUE 

MATIN077 


DO 100 1=1 *N 

MAT IN078 


DO 80 J=I,N 

MATIN079 


IF (IROW(J)-I) 80,85,80 

MAT IN080 

80 

CONTINUE 

MATIM081 

85 

IF (I-J) 90 » 100 *90 

MATIN082 

90 

DO 95 L=1 ♦ N 

MAT IN083 


LI=(I-1)»NDIM+L 

MAT IN084 


LJ=(J-1)*NDIM«-L 

MATIN085 


TEMP=A(LI) 

MATIN086 


A(LI>=A(LJ) 

MATIN087 

95 

A (L J) =TEMP 

MATIN088 


IROWI J)=IROWCI) 

MAT IN089 

100 

CONTINUE 

MATINO90 


DO 125 1=1, N 

MATIN091 


DO 105 J=I »N 

MATIN092 


IF (ICOL(J)-I) 105*110,105 

MAT IN093 

105 

CONTINUE 

MAT INQ94 

110 

IF (I-J) 115,125,115 

MAT IN095 

115 

DO 120 L=1,N 

MAT IN096 


IL=(L-1)*NDIM+I 

MAT IN097 


JL=(L-1 )*NUIM+J 

MATIN098 


TEMP=A(IL> 

MATIN099 


A(IL)=A(JL) 

MATIN100 

120 

A(JL)=TEMP 

MATIN101 


ICOL(J)=ICOL(I) 

MATIN102 

125 

CONTINUE 

MATIN1 03 


IROW ( NP1 ) =0 

MATIN104 


RETURN 

MATINl 05 


END 

MATIN106 
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C INPUT DATA 


TEST CASE 

10C SPHERICAL TANK 

3/4 FULL 

BOND 

NUMBER = 10 



12 12 

' 200 

0 1.8 

1.3 

.000001 

1.0 

10. 


0.0 

0.0 

.1216 

.0047 

.2438 

.0194 

.3665 

.0464 

.4874 

.0895 

.6008 

.1535 

.6960 

.2429 

.7319 

.2978 

.7564 

.3594 

.7657 

.4274 

.7527 

.5024 

.7209 

.5576 

.768 

.510 

.825 

.436 

.884 

.339 

.933 

.233 

.970 

.108 

.995 

-.037 

.999 

-.181 

.977 

-.341 

.925 

-.510 

.855 

-.647 





0.0 

-1.131 

.094 

-1.126 

.180 

-1.113 

.264 

-1.093 

.343 

-1.070 

.422 

-1.036 

.500 

-0.997 

.572 

-0.949 

.640 

-0.898 

.695 

-0.845 

.735 

-0.805 

.765 

-0.774 

0.0 

-.056 

0.0 

-.116 

0.0 

-.181 

0.0 

-.273 

0.0 

-.369 

0.0 

- .484 

0.0 

-.581 

0.0 

-.702 

0.0 

-.821 

0.0 

-.976 





1.0 

-65.788 

5.0 






-.36708 

3. 

14159267 

1.0 

0. 

0 
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INPUT DATA 


TEST - 

CASE 11 SPHEROIDAL TANK 

3/8 FULL 

BOND 

number = 5 



9 

12 200 

n 1,8 

1.3 

. OOOOQl 

i . r. 

5, 


0,0 

0,0 

.164 

,008 

.325 

, 032 

. 483 

.075 

,633 

,140 

, 768 

.231 

.827 

. 2ft7 

.8/7 

. 351 

,917 

,422 

.942 

,500 

.951 

.581 

,939 

.663 

,964 

.*00 

, 98 fi 

,54 0 

. 9941 

, 454 

1,00 

.367 

, 994 

,280 

, 98n 

. 194 

.954 

. Ill 



Q , 0 

- ,499 

,151 

-.489 

.301 

- . 460 

,447 

-.407 

.588 

-.333 

.718 

-.235 

.750 

- .201 

. 80 U 

-.154 

,833 

-,114 

.865 

-.070 

.900 

- , 015 

,923 

. 033 

0,0 

= ,06 

0.0 

-.12 

0 . 0 

-.18 

0 , 0 

- . 24 

o', 0 

- , 30 

0 , 0 

-.36 

n . o 

-.42 



1,0 

-16,12 

5,0 






= ,'394 

1. 

35 

1.0 

0. 

0 
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APPENDIX V. 

LISTING AND INPUT SAMPLE OF SSHAPE 
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PROGRAM SSHAPh ( INPUT,0urPUT,TAREi>sIWpUT,TAPt6=uUTPUT) 

PROJECT 02-l846-n2 
AUGUST 1969 

*LCW GRAVITY LIQUID-VAPOR I NTERF ACF SHAPES I N 
AX I SYMMET R I C CONT A INfcRS * 

CEC 6 6 n n FORTRAN 
D IRENS I Of' Y(10),F(10) 

DIMENSION XOR(20n0) , ZOR < 20 0 0 ) , YgX < 2 0 0 ) , Y BZ ( 2 0 0 ) 

DIMENSION T8DX(2),TBANG(2),A\S{3), ARG(3) 

D I PENS I ON XY8 ( 2 0 ft ) ,XRB(20O),XY(2 OOO),Xk(201)G) 

D IRENS ION' XKF(3),ZKF(3},XK8(3),ZKR(3) 

DI PENSION S(2onO),6(5o>,RS(5o>*SSw(5m,KN(2),RX(50J.Z(50) 
TYPE REAL Kfl,KUSV,KUP,K2F,K2B 
CONl = 1,74532925 e-02 
CON 2 = 57,2957795 
PI = 3,14159265 
NTAPEI = 5 
NTAPEO = 6 

C*****R E A D CONTAINER COORDINATES 
20ft0 READ (NTAPEI, 50fl) Np 
5f)0 FORMAT ( 2 1 5 ) 

READ (NTAPEI, 20ft) ( YBX ( I ) , I =1 , NP ) 

■2fl0 FORMAT ( 8F10.0 ) 

READ ( N!T APE 1 , 2ftft ) ( YBZ (I ) , I = 1 , NP ) 

C*****C .CNTAINER VOLUME 
NP1 = NP-1 
VT = 0, 

DO 95 I s 1 , N P 1 

VT s VT+,5*(YBX( I )**2 + Y8X( Y8Z( I )-YBZ( I+D ) 

95 CONTINUE 
VT = P I *vT 

2 0 READ (NTAPEI, 20ft) ALPHA , K 0 , Y 0 , BN , BETaD , DKG , DEC , DYQ 
READ (NTAPEI, 2 0 n > THET A , DTHET A , RLGTH 
C*****N ’UMBER OF EQUATIONS 
READ (NTAPEI, 50ft) Nl,jBOPT 
READ (NTAPEI, 500) NN(1),NN(2> 

DO 94 1=1, NP 
YBX(I) = Y8X( I >/RLGTH 
YBZ < I ) = YBZ ( I )/RLGTH 
94 CONTINUE 

Yft = YO/RLGTH 
DYO = DYO/RLGTH 
TPRINT = 2, 

I PR I NT =( (TPRINT + DTHETA/lO, )/ DTHET A) 

KftSV = NO 

DKCSV = DKO 
DKO = 4,*DKn 
YBX1 = Y8X ( 1 > 

YBZ1 = YBZ ( 1 ) 

II = 0 
JJ = 0 
JK = 0 

TBCX ( 2 ) = Kft 
10 no NT = o 

Y(l) = YO 
Y( 2) = 0, 

X = THET A*C0N1 
H s DTHET A*C0N1 
WRITE (NTAPEO, 3 n 0 ) 

3ft0 FORMAT (22H1 INPUT DATA) 

WRITE (NTAPEO, 3ft5) ALPH A , KO , Y 0 , BN 
3D5 FORMAT ( 9HOALPHA = ,ElS,8,6h ( DEG ) , 4x , 5HK0 = , El5 , 8, 6x , 5Hyo = 


SSH002QQ 
SSH00300 
SSH00400 
SSHQ0500 
SSH00600 
SSH00700 
SSH00800 
SSH00900 
SSH01000 
SSHOllO 0 
SSHO12Q0 
SSH01300 
SSHOHOO 
SSH01500 
SSH01600 
SSH01610 
SSH01620 
SSH01700 
SSH0180Q 
SSH01900 
SSH02200 
SSH02300 
SSH024Q0 
SSH0250 0 
SSH02600 
SSH027ft0 
SSH0280Q 
SSH02900 
SSH03000 
-SSH0310 0 
SSH03200 
SSH033O0 
SSH03400 
SSH03500 
SSH03600 
SSHG370Q 
SSH03800 
SSH03900 
SSH04000 
SSH04100 
SSH04200 
SSH04300 
SSH04400 
SSH04500 
SSH04600 
SSH04700 
SSHQ4800 
SSH04900 
SSH 05 0 0 0 
SSH05100 
SSH05200 
SSH05300 
SSH0540O 
SSH05500 
SSH0560D 
SSH05700 
SSH05800 
SSH05900 
SSH06000 
SSH0610 0 
SSH06200 



V-3 



1 El5 • 8 , 4 X , 5HBN = ,E15.8) 

SSH06300 


WRITE (NTAPEO,310) 

SSH06400 

3lo 

FORMAT (1HO,5X,5HTHETA,10X,3HXOR,12X J 3HZOR, 12X, 4 MY(1) ,12X, 

4hY(2)/)SSH06500 


XOR ( 1 ) = 0. 

SSH0660Q 


ZOR(l) s YO 

ssno67on 


XDEG = 0, 

SSh q68 p o 


WRITE! NTAPE0»315) XDEG # XOR ( 1 ) , ZCR ( 1 ) , Y < 1 ) * Y ( 2 ) 

SSH0690Q 


vui = 0 , 

SSH o7 0 0 0 


ICCUN7 = I 

SSH07100 


K = 1 

SSH07200 


IJK = 2 

SSHo73qo 


00 35 1=2, nr 

SSH074Q 0 

75 

F C 1 ) = Y { 2 ) 

SSH07500 


IF ! X , EQ , 0 , ) 125.130 

SSH07600 

125 

F(S) s Y 0 * ( 1 , - K 0 > 

SSH077 oo 


GO TO 135 

SSH07800 

130 

CONTINUE 

SSH07900 


TERM = Ytl>**2+Y<2)**2 

SSH08000 


IF < TERN-1, QE+6Q ) 131,132,132 

SSH 081 0 0 

132 

XDEG s X*C0N2 

SSH082QG 


WRITE ( NT AP£0 , 3l5 ) XDEG, Y < 1) , Y ( 2) , F ( 1 ) , F (2) 

SSH083Q0 


IF { JJ,EO.O ) 777,778 

SSH08400 

777 

KO s Ko-DKOSV 

SSHq 85 o o 


GO TO 1000 

SSHO860O 

778 

CONTINUE 

SSH08700 


GO TO 74 

SSH08800 

131 

CONTINUE 

SSH 089 0 0 


F ( 2 ) = ( C2.*YC1)**2 + 3,*Y ( 2 ) ** 2 ) / Y ( 1 ) )-( Y (2>/Y (1)**2)*(1,/ 

SSH0900Q 


1 TAN(X) )*(Y(l)**2+V(2)**2) + (l,/Y(l))*(BN*(yU) 

2 <{2,*Kn)/Yn>>*<SQRT<Y<l)**2+Y<2)**2) )**3 
135 CONTINUE 

S = RKLDEQ ( N, Y, F, X, H, N T ) 

IF < S-1,0 ) 105,75,110 


★COb(X)-YO)- SSN09100 
SSH09200 
SSH0930 0 
SSH09400 
SSH09500 


105 STOP 

110 TERM1 = SIN(X) 

TERM2 = COS(X) 

XDEG = X* CON2 

TERM = YBZ< I -»! ) -YB2 ( I ) 

IF { TERM.EO.n, )465,470 
465 R = 1.0E+10O 
GO TO 475 

47 o R = (YBX< I)-YBX( 1-1) )/TERN 


SSM09600 

SSH097Q0 

SSH09800 

SSHG9900 

SSH10000 

SSH10100 

SSH10200 

SSH10300 

SSH10400 


475 CONTINUE SSH10500 

YX = <YBX( I-1)+R*Y8Z< 1-1) )/ CTERM1+R*TERM2) SSH10600 

YBX2 = YBX ( I ) SSH1070Q 

YBZ 2 = YB2 { I > SSH1Q600 

YZ = SQRT(YBX( IJ**2+YdZ( I)**2) SSH10900 

K s K+l SSHllOOO 


XOR { K ) a Y(1)*TERM1 
20R ( K ) s Y(D*TERM2 

VU1 = VUl+,5*(XOR(K-l)*»2+XOfi(K)**2)*(ZOR(K-l)-ZOR<K))*Pl 
IF ( YZsfeO.O, ) 111,113 

111 THETAB s 90,*CON1 
GO TO 112 

113 THETA8 = AS I N ( Y8X ( I ) / YZ ) 

112 CONTINUE 

IF < Y ( 1 ) * YX ) 40,45,45 
40 CONTINUE 

IF ( XDEG«=90, ) 50,74,74 
50 CONTINUE 

IF ( I COUNT *■ I PR I NT ) 115,120,115 
115 ICOUNT = ICOUNT+1 


SSH11100 
SSN11200 
SSH113Q0 
SSH11400 
SSH11500 
SSH116 0 0 
SSH1170Q 
SSHH 800 
SSH11900 
SSH12000 
SSH121Q0 
SSH12200 
SSH12300 
SSH124 0 0 
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GO TO 117 

SSH12500 

120 

CONTINUE 

SSH126 0 0 


WRITE { NTAPEO, 3i5 ) XDEG, XOR ( K ) , ZOK ( K) , y < 1 ) , Y ( ? ) 

SSH12800 

315 

FORMAT C 0E15.7 ) 

SSH12900 


1CCUNT = 1 

SSH13000 

117 

IF ( X-ThETAB ) 116,118,118 

SSH13200 

116 

YBX1 = YBX2 

SSH1330Q 


YBZl = YBZ2 

SSH134Q0 


GO TO 75 

SSH13500 

11 a 

ISAVE = I 

SSH1360 0 


GO TO 35 

SSH137 0 0 

45 

GAMMA =(ATAN2(Y(2)*SIN(X)+Y(1)*00S(X),Y( 1 ) *S I N ( X ) - Y ( 2 ) * COS ( X > > )* 

SSH13800 


1 CCN2 

SSH13900 


IF < ABS ( GAMMA ) -9n , ) 46,46,74 

SSH14000 

46 

CONTINUE 

SSH14100 


ARG1 = -(YBX(I )-YBX< 1-1) ) 

SSH14200 


ARG2 = - (YBZ ( I )-Y8Z(I-l)) 

SSH143Q0 


IF ( ABS(ARG2) ,LE t (1.0E-08) ) 47,48 

SSH14400 

47 

PHI s 90, 

SSH14500 


GO TO 49 

SSH14600 

40 

PHI = ATAN2{ARG1,ARG2)*C0N2 

SSH14700 

49 

ANGLE = GAMMA+PHI 

SSH14800 


WRITE ( NT APEO, 3l5 ) XDEG , X OR ( K ) , ZOR ( K ) , Y ( 1 > , Y ( 2 ) 

SSH14900 


WRITE ( NT APEO, 3 2 0 5 GAMMA, PHI , ANGLE, Yd, KQ 

SSH150G0 

320 

FORMAT (IQHfl GAMMA = ,E15.8,8H PHI = ,El5.8,10H ALPHA = ,E15,8, 

SSH15100 


1 7'H Y0 = , £15 , 8, 7H K0 = ,E15,8> 

SSH152 0 0 


IF ( ABS< ANGLE-ALPHA)-, 5 ) 55,55,6n 

SSH153 0 0 

55 

VU2 = 0, 

SSH15400 


JJ = 0 

SSH155 0 0 


JK c 0 

SSH15600 


TBDXI2) = KASV 

SSH1570 0 


YBXSV = YBX I JSAVE-1) 

SSH15800 


YB2SV s YBZ ( ISAVE-1 ) 

SSH1590Q 


YBX(ISAVE-l) = XOR { K ) 

SSH16000 


Y B 2 ( ISAVE- 1} = ZOR(K) 

SSH16100 


IJ = ISAVE-1 

SSH162 0 0 


DO 100 I J, NP1 

SSH16300 


VU2 s VU2+,5*(YBX(J)**2 + Y8X(vI+1)**2)*{YBZ( J)-YBZ( J+1))*PI 

SSH164 0 0 

ino 

CONTINUE 

SSH165 0 0 


ybx n Save* 1 ) = ybxsv 

SSH16600 


YBZ ( ISAVE-1) = YBZSV 

SSH167 0 0 


VUl = VU1*RLGTH**3 

SSH168QQ 


VU2 = VU2*RIGTH**3 

SSH16900 


VU = VU1+VU2 

SSH17000 


VL c VT»VU 

SSH171QQ 


BETA s vu/vt 

SSH17200 


WRITE < NTAPEO, 325 ) VUl, VU 2, VL , VL , V T, BET A 

SSH17300 

i325 

'FORMAT (8H0 VUl = ,El5,8,8H VU2 = ,E15,8,7H VU = #El5 t B» 

SSH17400 


1 7>H VL s , E15 ,8, 7H VT = ,£15.8y/9H BETA = ,E15,8j 

SSH17500 


BB = BETAD-BETA 

SSH176Q0 


IF ( IBOPT.EQ.O ) GO TO 145 

SSH17700 


IF ( A8S(BB)*D8C ) 145,145,140 

SSH17800 

,14q 

CONTINUE 

SSH17900 


IF ( II ) 155,150,155 

SSH18000 

15 0 

II s II+l 

SSH181QQ 


ANS(l) 5 Y0 

SSH18200 


YO * YO+DYO 

SSH18300 


ARGU) s BETA 

SSH104OO 


K n = KOSV 

SSH18500 


DK0 = DKOSV 

SSH18600 


GO TO 1000 

SSH1B700 

155 

IF ( 11*1 ) 165,160,165 

SSH18800 
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160 II = II + l 

ANS C 2 ) = YO 
ARC C 2 > = BETA 
Y 0 = YO + DYO 
K 0 a KOSV 
DKO = DKOSV 
GO TO 1000 

165 IF ( I I * 2 ) 175,170, 175 
170 ANS { 3 ) s YO 
II = II+l 
175 ARG { 3 ) s BETA 

IF < ARG ( 1 ) - ARG ( 3 ) ) 405,400,400 
400 AT = ARG ( 1 ) 

ARG(l) = ARG ( 3 ) 

ARG ( 3 ) 3 AT 
7S * ANS(l) 

.ANS ( 1 ) = ANS ( 3 ) 

ANS < 3 > a TS 

405 IF < ARG ( 1 ) • ARG ( 2 ) ) 415,410,410 
4lo PP = ARG(l) 

ARG(l) = ARG(2) 

ARG ( 2 ) = PP 
TS * ANS ( 1 ) 

ANS(l) s ANS<2) 

ANS ( 2 ) = TS 

415 IF < ARG(2)wARG(3> ) 425,420,420 
42o -ST = ARG(2) 

ARG < 2 ) = ARGO) 

ARG ( 3 ) s ST 
TS s ANS ( 2 } 

ANS < 2 ) = ANS 1 3 ) 

ANS ( 3 ) s TS 
•425 .CONTINUE 

IF < 8ETAD-ARGC1) > 435,430,430 
43 o IT < BET AD- ARG C 3 ) ) 44 0 ,445,445 
44o IF < BET AD- ARG ( 2 ) ) 450,455,455 

450 -ANSE = ( ANS ( X ) * ( ARG ( 2 > -BETA!) ) -ANS < 2 ) * ( ARG ( 1 ) -BET AD) >/(ARGl2)- 
1 ARG(l)) 

GO TO 46o 

455 ANSE = (ANS(2)*(ARG(3>’*BETAD)-ANS(3)*(ARG(2)’ , 8ETAD))7(ARG(3)- 
1 ARG ( 2 > ) 

GO TO 46 o 

435 ANSE 3 (ANS(1)*(ARG(2)’*BETAD)+ANS(2)*(SETAD-ARG(1) > )/(ARG(2)- 
1 ARG(l)) 

GO TO 46o 

445 ANSE = (ANS(2)*CARG(3)-BETAD)+ANS(3>*C8ETAD*ARG<2> >)/(ARG(3)«- 
1 ARG< 2) ) 

460 'CONTINUE 

WRITE (NTAPE0,7H0) ANSE 
700 FORMAT (2QH0 EXTRAPOLATED YO s ,615,8) 

YO e ANSE 
ARG ID = ARG < 2 ) 

ARG { 2 ) s ARG { 3 ) 

ANS ( 1 ) s ANS ( 2 ) 

ANS ( 2 ) s ANS ( 3 ) 

ANS ( 3 ) s YO 
DKO = DKOSV 
Kfi s KOSV 
GO TO 1000 
145 CONTINUE 
SCI) = 0, 

•SSV(l) s SCI) 


SSH1890 0 
SSH1900Q 
SSH19100 
SSH19200 
SSH19300 
SSH194 0 D 
SSH1950 0 
SSH196 0 0 
SSH19700 
SSH19900 
SSH19900 
SSH20000 
SSH20100 
SSH2 o2 0 0 
SSH2q300 
SSH20400 
SSH2q500 
SSH2q600 
SSM20700 
SSH2080Q 
SSH2Q900 
SSH210 0 0 
SSH21100 
SSH21200 
SSH21300 
SSH21400 
SSH21500 
SSH21600 
SSH21700 
SSH218O0 
SSH219.00 
SSH22000 
SSH22100 
SSH22200 
SSH223Q0 
SSM22400 
SSH22500 
SSH22600 
SSH22700 
SSH2280 0 
SSH22900 
SSH23000 
SSH23100 
SSH232Q0 
SSH233 0 0 
SSM23400 
SSH235qo 
SSH23600 
SSH23700 
' SSH230OO 
SSH239Q0 
SSH24000 
SSH24100 
SSH242Q0 
SSH24300 
SSH2440Q 
SSH245Q0 
SSH2460 0 

SSH24700 

SSH24800 

SSH24900 

SSH25000 
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NS = K-l 
DO 1 L*2,MS 

DEUX s <X0R(U+1)-X0R(L) )*RLGTh 
DELY = I20RUW0R(L + 1> )*RU»Tl- 
CEL5 s SGRT( DELX**2 + t}ELY**2) 

S(U S S(L-1) + DELS 

1 continue 

DO 5 J [ si , 2 

IF INN-(JI) , fcG.G ) GO TO 6 
NN1 = NN(JI)-1 
NM s NNI(UI) 

SI s S(NS)/NtMl 
G ( 1 ) = 0, 

Rsa> * i, 

DO 2 Ls 2,NN1 
DO 3 J s l> NS 

IF (SI J)-SI*(L-1) ) 3,4,4 

4 SSV ( L ) = SIJ) 

HI ? S(Jfl)-SIJ) 

H2 = S(J*1)-S(J) 

DEN = Hl*H2* < H2-H1 ) 

APO = -IH2**2-Hl**2)/D£N 
API = h 2**2/D£N 
APS = * I Hi**2 ) /DEN 
APPO = 2 , / ( H1*H2 ) 

APF1 = ( *2 , * H2 ) / DEN 
APr2 = (2,*h1)/DEN 
RXIL) = XOft ( J ) *RLGTH 
Z II) = (YO-ZORI J))*RLGTH 

RPS =( AP2*X0RI J+l)+AP0 *XOR( J)+aP 1*XOR( J-l) )*PL6TH 

RPSS =(APP2*X0R( J + 1)+APP0*XUR( J)+APP1*X0R( J-l) )+RI_GTH 

ZPS =( AP2*(YQ-ZQR( J+i) ) + AP 0* I Y Q'-ZOR I J5 ) + API* I Y 0-ZQR U-l ) 3 )*RLGTH 

ZRSS =( APP2*{ YO-ZORI J+U ) + APP 0+ I YO-ZGR < J) ) ♦ APPl* ( Y 0-ZQR ( J- 1 ) ) >* ' 

1 rlgtk 

CKAP2 s zps/ ixori J)* ftLGTH> 

CKAP1 s RPS*ZPSS-ZPS*RPSS 
GIL) = CKAP2**2+CKAPi**2 
RS(t) = RPS 
GO TO 2 
.3 CONTINUE 

2 .CONTINUE 
SSVINH) s SINS) 

GI-NM) 5 G(NNl) 

RS(NM) = RS ( NNl } 

RXINM) = XOR(K)*RLGTH 
ZINM) = (Y0-2ORIK) )*RLGTH 
PRINT 314 

314 FORMAT{lHl,l4x,*G*,2oX,*RS*,19X##S*,l9X,*XR*,l9X i *XY*/) 

WRITE ( NT APEO » 319 ) II , G < 5) , HE 1 1 ) , SSV (I) . RX I I) , Z 1 1 ) , 1 = 1. ) 

319 FORMAT I l5,5E2n,8 ) 

5 CONTINUE 

6 CONTINUE 

YBX(iSAVE-l) = Y8XSV 
YBZ I I SAVE” 1 ) * Y8ZSV 
DO 605 «J = 1 t NP 
XRBIJ) s YBXIJ)*RLGTH 
XYEU) siYO-VBZI J) )*RLGTH 
605 CONTINUE 

XKF{1) = XOR I K* 2 ) 

XKF I 2) = XOR(K-l) 

XKFI3) s XORIK) 

ZKFU) a f-ZOR I K = 2 ) +ZOR 1 1 ) 


SSH25100 
SSH25200 
SSH253Q0 
SSW25400 
SSH25500 
SSH2560Q 
SSH25700 
SSH2580Q 
SSH25900 
SSH26000 
SSH261QQ 
SSH26200 
SSH26300 
SSH264 o 0 

SSH26500 
SSH26600 
SSH267 0 0 
SSH268 0 0 
SSH26900 
SSH27000 
SSH27100 
SSH2720 0 
SSK273oo 
SSH274QQ 
SSH27500 
SSH276Q0 
SSH27700 
SSH27800 
SSH279Q0 
SSH28000 
SSH2810 0 
SSH28200 
SSH28300 
SSH2840 0 
SSH28500 
SSH28600 
SSH28700 
SSH28800 
SSH2S900 
SSH290 0 0 
SSH291GQ 

SSH29200 
SSH293 00 
SSH294D0 
SSH295Q0 
SSH2960D 
SSH29700 
SSH29B00 
SSH29900 
SSH30000 
SSH3Q100 
SSH30200 
SSH303OO 
SSH30400 
SSH30500 
SSH3 06 0 0 
SSH30700 
SSH30BOG 
SSH3()9O0 
SSH31OO0 
SSH31100 
SSH3i200 
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ZKF<2) =-20R(K-l )+ZOR(l ) 

ZKF<3) =-ZOR( K)+ZOR( 1) 

CALI SFIT ( YR,XK3,XK, YR2.XKI ,Z«F ) 

K2F = A8S(XK ) 

FHR2 = K2F*(1 i +YR2**2>**(1,5) 

XKB(l) = Y8X ( I -1 ) 

XKB ( 2 ) = XORCK) 

XKB ( 3 ) = YBX( I ) 

ZK8(1) YBZU-D+ZORU) 

ZKR ( 2 ) ZOR(K)+ZOR(l) 

ZKE(3) s- YBZ < I ) + ZOR ( 1 ) 

CALI SFIT <YR,XK3,XK,YR2,XKB,ZKB) 

K2B = ABSCXK > 

ZRR2 s k2B*(1,+YR2*+2}**<1,5> 

CGAM = (1 ,/SlN(ANGLE*CONl) )*( K2B* COS C ANGLE*COM ) * K2F) 

WRITE (NTAPEO,Bn5) «2F, FRR2 , k2r, Z«r2, CGA* 

8n5 FORMAT <9H0 K2F = ,El5,8,9H FRR2 = ,6l 5 .S//7H K2B = *1=15,8, 
1 9F ZRR2 = .E15.8//8H OGAM = ,Ei5,8) 

DO 606 Jal,K 
XR ( J)= XOR ( v) *RLGTH 
X Y ( J)=(YO“ZOR( J) )*RLGTH 
606 CONTINUE 

YO = Y0*RLGTH 
WRITE <NTAPE0,316> YO 
318 FORMAT (7H0 Yn = ,El5,8) 

WRITE ( NT APEO, 335 ) 


SSH31300 
SSH3l 4 g 0 

SSH315Q0 
SSH31600 
SSH31700 
SSM31800 
SSH31900 
SSH32000 
SSH321Q0 
SSH32200 
SSH32300 
SSH3240O 
SSH32500 
SSH32600 
SSH327 0 0 
SSH328Q0 
SSH32900 
SSH33000 
SSH33100 
SSH33200 
SSH3330 0 
SSH33400 
SSH33S00 

SSH33600 

SSH33700 

SSH33QQ0 


335 FORMAT (1 HI , 9X , 3HXRB , 17X , 3riX Y 6/ ) SSH33900 

WRITE INTAPE0,317) ( XRB ( J ) , X YB ( J ) * J= 3 , N'P > SSH34qoo 

317 FORMAT { 2£20,8 > SSH34100 

SUM 1 = 0, SSH34200 

DO 610 J=2,K SSH34300 


SUMl = SUMl+,5*<XR( J)**2*XY( J ) +XR (J- 1 ) ** 2*x Y < J-l) )*(XY( J>*XY( J-l) )SSH344 qo 


610 CONTINUE 
SUM2 =■ 0, 
IJ = ISAVE 


SSH345QG 

SSH34600 

SSH34700 


XRBSV s XRB < I J) 

XYESV s XYB(IJ) 

XRB ( I J ) = XOR(K)*RlGTH 
XYBUJ) = (YO-ZOR(K) )*RLGTH 
DO 615 J= I J » NP1 

SUM 2 = SUM2+,5*( XRB< J)**2*XYB( J)+XR8< J+1)**2*XY8< J*U )*CXYB( J+l)- 
1 XYB< J) ) 

615 CONTINUE 

BZCGU s ( RI*(SUM1 + SUm2) )/VT 
XYB(IJ) = XYBSV 
XRB < I J > s XRBSV 
SUMl = 0, 

DO 620 J=1,NP1 

SUM.l = SUM1+,5*CXR6{ J)**2*XYB( J ) + X RB ( J+1)**2*XYB< J + l) )*(XYB( J+l) - 
1 XY8 ( J) ) *P I 
620 CONTINUE 

ZCGT s SUM1/VT 

ZCG = <ZCGT-8ZCGU)/(1,-BETA) 

WRITE (NTAPtO, 8nQ) ZCGT , BZCGl , ZCG 

800 FORMAT ( 9H0 ZCGT = ,El5,8,10H BZCGU = ,£15,8,88 2CG = ,El5,8) 

GO TO 2000 
60 CONTINUE 

IF ( ANGLE-ALPHA ) 74,65,71 

71 IF ( JJ*1 ) 72,73,72 

72 TBDX(l) = KO 
TBANG(l) = ANGLE 
JJ = JJ + 1 


SSH3480D 
SSH34 9 o o 
SSH35000 
SSH35100 
SSH352 0 0 
SSH35300 
SSH35400 
SSH35500 
SSH35600 
SSH35700 
SSH35800 
SSH35900 
SSH360 0 0 
SSH36100 
SSH36200 
SSH36300 
SSH364Q0 
SSH36500 
SSH366GQ 
SSH367 go 
SSH36800 
SSH36900 
SSH37QO0 
SSH371Q0 
SSH37200 
SSH373 0 0 
SSH3740Q 
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GO TO 70 

SSH37500 

73 

TBEX<2) = ko 

SSH376 o 0 


TB ANG< 2 ) = ANGLES 

SSH37700 


IF < JK > 76,76,77 

SSH378O0 

77 

K n = , 5* { K QP+TBDX < 2 ) ) 

SSH379 Q 0 


GO TO 1000 

SSH38000 

76 

KO = (TBDX(l)*(T8ANG<2>’-ALPHA)-T8DX(2)*(TBAN&(l)-ALPHAn/ 

SSN38100 


1 ( TBAIMG(2)-T6ANG(1> ) 

SSH382Q0 


JJ = 1 

SSH3B3 o 0 


TBCX(l) s T8DX ( 2 ) 

SSH38400 


TBANG(l) = T8ANG { 2 } 

SSH38500 


GO TO 1000 

SSH386Q0 

65 

DKO = DKO/4, 

SSH38700 


KO s KO - DKO 

SSH388Q0 


GO TO 1000 

SSH38900 

70 

DKO = DK0/4, 

SSH39000 


KO = Kq + DKQ 

SSB39100 


GO TO 1000 

SSH392Q0 

74 

KOf = KQ 

SSH3930 0 


Kn s , 5* { K0 + T8DX ( 2 ) ) 

SSH39<00 


JK s 1 

SSH395Q0 


GO TO 1000 

SSH396Q0 

35 

continue 

SSH39700 


end 

SSH3980 0 
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FUNCTION RKLDEQ(N,Y,F,X,H,N1 ) RKLDQ T 

D2 UCSC RKLUEG RUNGE-KUTTA-GILL LINEAR DIFFERENTIAL FQUaTION bOLVER F-63 X 
C2 UCSO RKLDEQ F 63 

MODIFIED NAY 1963 <Q REMOVED FRUM CALLING SfcLUtNCE) 


TEST OF ALGOL ALGORITHM 


c 

DIMENSION Y(ln),F(lQ),G(lO) 

REAL X, H-- INTEGER N, NT- -COMMENT-- BEG I N INTEGER I,J,L-REAL A 

RKLDQ 

2 


N T = N T + 1 

RKLDQ 

3 

c 

GO TO (1,2, 3, A ), NT 
GO TO S ( NT ) 

RKLDQ 

4 

1 

DO 11 J=1,N 

RKLDQ 

5 

11 

O('v') = 0, 

RKLDQ 

6 


A=,5 

RKLDQ 

7 


XsX+H/2, 

RKLDQ 

8 


GO TO 5 

RKLDQ 

9 

2 

A=, 29289321881 

RKLDQ 

10 


GO TO 5 

RKLDQ 

11 

3 

A = 1 ,7071067812 

RKLDQ 

12 


X=X+H/2, 
GO TO 5 

RKLDQ 

13 


RKLDQ 

14 

4 

DO 41 I=1,N 

RKLDQ 

15 

41 

YU ) = Y ( I)+H*F( I)/6,*-0( I )/3, 

RKLDQ 

16 


NT= 0 

RKLDQ 

17 


RKLDEQs2, 

RKLDQ 

18 


GO TO 6 

RKLDQ 

19 

5 

CO 51 L=1.N 

RKLDQ 

20 


Y(L)sY(L)+A*<h*F(L)-Q{L)) 

RKLDQ 

21 

51 

0{L)*2,*A*H*F(L)+(1,-3.*A)*Q(L) 

RKLDQ 

22 


RKLDEOsl, 

RKLDQ 

23 

6 

CONTINUE 

RKLDQ 

24 


RETURN 

RKLDQ 

25 


end 

RKLDQ 

26 
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SUBROUTINE SFIT<YR,XK3,XK,YR2,R,YF> 

FIT QUADRATIC THROUGH THREE POINTS 
DIMENSION R ( 3 ) , Y T ( 3 ) 

IF(ABSt (YTC3)-YT<2) )/(R(3>-R(2) ) ) ,GT.1.0) CO TO 30 
IF(ABS< (YT(2)-YT (1) )/(R(2)-H(l!) ) ,GT.l, 0) GO TO 3o 

XTEMP=(R{1)**2-K(3)**2)* (R(2)-R(3))-(R(2)**2-RC3)**2)*(RU)»R(3)) 
YTEMP=( YT ( 1)--YT(3 ) )*<R(2)-R(3n-(YT(2)-YT(3n*(R<l)-R(3)) 
ADsYTEMP/XTEHP 

AB=(YT(2)-YT(3)-AU*(R(2)**2-R(3)**25)/(R(2)-R(3)) 

YR=2 f 0*AD*R(3)+A3 
YR2=2,o*AD*R<2)+AB 
YRR=2, 0*AD 

XK=YRR/CL, 0 + YR2**2>**l,5 
H2 sR < 2 ) fjR ( 3 ) 

HlsR(l)-R( 3) 

F2=YT(2> 

Fls YT ( 1 ) 

Ffl = YT ( 3 ) 

YRRa2,*{H2*(Fl-F0)-Hi*(F2-FO) ) / < H1+H2* < HI -H2 > ) 

XK35YRR/(i.n+YR*+2)**l,5 

RETURN 

3 0 YTEMP={ YT < 1)**2-YT(3)**2)*<YT(2)..YT(3>)-(YT(2)**2-YT{3)**2)* 

1 ( YT < 1)*YT(3) ) 

XTEMPs{R(l)-H(3) )*(YT(2)-YT(3) >-<R(2)-R(3) )*CYT(J )“YT(3M 
ADsXTEMP/YTEMP 

AB= (R(2)-R(3)-AD* (YT<2)**2-YT(3)**2) >/<YT(2)-YT(3)) 

YR*1, 0/(2, 0*AD*YT<3)+AB) 

YR2=l,0/(2.n*AD*YTt2)+AB) 

YRR=-2 , 0*AD*YR2**3 
XKsYRR/(l, 0+YR2**2>**1.5 
H2 = YT < 2 ) - YT < 3 j 

H1=YT(1)-YT<3) 

F2=R<2) 

F1 = R < 1 ) 

Fn = R { 3 > 

RYY=2,+<H2*<F1-F0)-H1*<F2-F0) ) / ( H1*H2* C H1-H2 ) ) 

YRR3=-RYY*YR**3 
XK3sYRR 3/( 1, 0+YR**2)**l,5 

RETURN 
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